1) This must be a dense matrix, so you could directly copy in the values with a strided loop<div><br></div><div>2) Is there a reason for the output to be a matrix, rather than a collection of vectors?</div><div><br></div><div>
   Matt<br><br><div class="gmail_quote">On Tue, May 18, 2010 at 2:52 AM, Umut Tabak <span dir="ltr">&lt;<a href="mailto:u.tabak@tudelft.nl">u.tabak@tudelft.nl</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Dear all,<br>
<br>
I would like to do a Gram-Schmidt orthogonalization of a matrix with say 10 columns, I wrote a function for this, However, I do not do this in place for the moment. I create a new matrix for the orthogonal(or M-orthogonal columns). And set the columns of that new matrix, I could not see how I could do this in place, where I read the original matrix in binary format from a file. Then its assembly is finalized I suppose, is there still a way to make changes on that matrix in place?<br>

<br>
The code is below.<br>
<br>
A simple translation of a matlab script(more efficient implementations are possible I am sure.)<br>
<br>
int mgrasch_Morth(Mat &amp;inMat, Mat &amp;outMat, Mat M )<br>
{<br>
 // read size of the matrix<br>
 PetscErrorCode ierr;<br>
 PetscInt r,c;<br>
 ierr = MatGetSize(inMat, &amp;r, &amp;c);                                          CHKERRQ(ierr);<br>
 ierr = MatDuplicate(inMat, MAT_COPY_VALUES, &amp;outMat); CHKERRQ(ierr);<br>
 // set up the necessary workspace<br>
 Vec Q,R,dummy1,dummy2,dummy3,orthoVec;<br>
 PetscScalar val1, val2;<br>
 //<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;Q);     CHKERRQ(ierr);<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;R);     CHKERRQ(ierr);<br>
<br>
 // dummy vectors for matrix multiplication<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;dummy1);     CHKERRQ(ierr);<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;dummy2);     CHKERRQ(ierr);<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;dummy2);     CHKERRQ(ierr);<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;dummy3);     CHKERRQ(ierr);<br>
 //<br>
 ierr = VecCreate(MPI_COMM_SELF, &amp;orthoVec);     CHKERRQ(ierr);<br>
 //<br>
 ierr = VecSetSizes(Q, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 ierr = VecSetSizes(R, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 ierr = VecSetSizes(dummy1, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 ierr = VecSetSizes(dummy2, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 ierr = VecSetSizes(dummy3, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 ierr = VecSetSizes(orthoVec, r, PETSC_DECIDE);  CHKERRQ(ierr);<br>
 //<br>
 ierr = VecSetFromOptions(R);         CHKERRQ(ierr);<br>
 ierr = VecSetFromOptions(Q);         CHKERRQ(ierr);<br>
 ierr = VecSetFromOptions(dummy1);    CHKERRQ(ierr); // for M * A(:,m)<br>
 ierr = VecSetFromOptions(dummy2);    CHKERRQ(ierr); // for M * Q<br>
 ierr = VecSetFromOptions(dummy3);    CHKERRQ(ierr); // for A(:,m) above<br>
 ierr = VecSetFromOptions(orthoVec);  CHKERRQ(ierr); //  orthogonalized result vec<br>
 // create the arrays for MatSetValues<br>
 // row<br>
 PetscInt rowArr[r];<br>
 for(int k=0;k!=r;k++)<br>
   rowArr[k]=k;<br>
 // to set one column at a time<br>
 PetscInt C = 1;<br>
 PetscInt colArr[1];<br>
 PetscScalar valsOrthoVec[r];<br>
 //<br>
 for(int i=0;i!=c;i++){<br>
   MatGetColumnVector(inMat, Q, i);<br>
   for(int j=i+1;j!=c;j++){<br>
     ierr = MatGetColumnVector(inMat, dummy3, j); CHKERRQ(ierr);<br>
    ierr = MatMult(M,dummy3,dummy1);             CHKERRQ(ierr);<br>
     // val1 -&gt; Q&#39;* M * A(:,m)<br>
     ierr = VecTDot(Q, dummy1,&amp;val1);             CHKERRQ(ierr);<br>
     ierr = MatMult(M,Q,dummy2);                  CHKERRQ(ierr);<br>
     // val2 -&gt; Q&#39;* M * Q<br>
     ierr = VecTDot(Q,dummy2,&amp;val2);              CHKERRQ(ierr);<br>
     //<br>
     ierr = VecAXPY(dummy3, -val1/val2, Q);       CHKERRQ(ierr);<br>
     ierr = VecCopy(dummy3, orthoVec);            CHKERRQ(ierr);<br>
     // set the orthogonal vector<br>
     colArr[0] = j;<br>
     //<br>
     ierr = VecGetValues(orthoVec, r, rowArr, valsOrthoVec);                                           CHKERRQ(ierr);<br>
     ierr = MatSetValues(outMat, r, rowArr, C, colArr, valsOrthoVec, INSERT_VALUES);  CHKERRQ(ierr);<br>
   }<br>
 }<br>
 ierr = MatAssemblyBegin(outMat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
 ierr = MatAssemblyEnd(outMat, MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);<br>
 return 0;<br>
}<br>
<br>
</blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>

</div>