[petsc-users] Advice on how to tackle with in place operation for a Matrix

[Matthew Knepley]

1) This must be a dense matrix, so you could directly copy in the values
with a strided loop

2) Is there a reason for the output to be a matrix, rather than a collection
of vectors?

   Matt

On Tue, May 18, 2010 at 2:52 AM, Umut Tabak <u.tabak at tudelft.nl> wrote:

> Dear all,
>
> 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?
>
> The code is below.
>
> A simple translation of a matlab script(more efficient implementations are
> possible I am sure.)
>
> int mgrasch_Morth(Mat &inMat, Mat &outMat, Mat M )
> {
>  // read size of the matrix
>  PetscErrorCode ierr;
>  PetscInt r,c;
>  ierr = MatGetSize(inMat, &r, &c);
>  CHKERRQ(ierr);
>  ierr = MatDuplicate(inMat, MAT_COPY_VALUES, &outMat); CHKERRQ(ierr);
>  // set up the necessary workspace
>  Vec Q,R,dummy1,dummy2,dummy3,orthoVec;
>  PetscScalar val1, val2;
>  //
>  ierr = VecCreate(MPI_COMM_SELF, &Q);     CHKERRQ(ierr);
>  ierr = VecCreate(MPI_COMM_SELF, &R);     CHKERRQ(ierr);
>
>  // dummy vectors for matrix multiplication
>  ierr = VecCreate(MPI_COMM_SELF, &dummy1);     CHKERRQ(ierr);
>  ierr = VecCreate(MPI_COMM_SELF, &dummy2);     CHKERRQ(ierr);
>  ierr = VecCreate(MPI_COMM_SELF, &dummy2);     CHKERRQ(ierr);
>  ierr = VecCreate(MPI_COMM_SELF, &dummy3);     CHKERRQ(ierr);
>  //
>  ierr = VecCreate(MPI_COMM_SELF, &orthoVec);     CHKERRQ(ierr);
>  //
>  ierr = VecSetSizes(Q, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  ierr = VecSetSizes(R, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  ierr = VecSetSizes(dummy1, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  ierr = VecSetSizes(dummy2, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  ierr = VecSetSizes(dummy3, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  ierr = VecSetSizes(orthoVec, r, PETSC_DECIDE);  CHKERRQ(ierr);
>  //
>  ierr = VecSetFromOptions(R);         CHKERRQ(ierr);
>  ierr = VecSetFromOptions(Q);         CHKERRQ(ierr);
>  ierr = VecSetFromOptions(dummy1);    CHKERRQ(ierr); // for M * A(:,m)
>  ierr = VecSetFromOptions(dummy2);    CHKERRQ(ierr); // for M * Q
>  ierr = VecSetFromOptions(dummy3);    CHKERRQ(ierr); // for A(:,m) above
>  ierr = VecSetFromOptions(orthoVec);  CHKERRQ(ierr); //  orthogonalized
> result vec
>  // create the arrays for MatSetValues
>  // row
>  PetscInt rowArr[r];
>  for(int k=0;k!=r;k++)
>   rowArr[k]=k;
>  // to set one column at a time
>  PetscInt C = 1;
>  PetscInt colArr[1];
>  PetscScalar valsOrthoVec[r];
>  //
>  for(int i=0;i!=c;i++){
>   MatGetColumnVector(inMat, Q, i);
>   for(int j=i+1;j!=c;j++){
>     ierr = MatGetColumnVector(inMat, dummy3, j); CHKERRQ(ierr);
>    ierr = MatMult(M,dummy3,dummy1);             CHKERRQ(ierr);
>     // val1 -> Q'* M * A(:,m)
>     ierr = VecTDot(Q, dummy1,&val1);             CHKERRQ(ierr);
>     ierr = MatMult(M,Q,dummy2);                  CHKERRQ(ierr);
>     // val2 -> Q'* M * Q
>     ierr = VecTDot(Q,dummy2,&val2);              CHKERRQ(ierr);
>     //
>     ierr = VecAXPY(dummy3, -val1/val2, Q);       CHKERRQ(ierr);
>     ierr = VecCopy(dummy3, orthoVec);            CHKERRQ(ierr);
>     // set the orthogonal vector
>     colArr[0] = j;
>     //
>     ierr = VecGetValues(orthoVec, r, rowArr, valsOrthoVec);
>                           CHKERRQ(ierr);
>     ierr = MatSetValues(outMat, r, rowArr, C, colArr, valsOrthoVec,
> INSERT_VALUES);  CHKERRQ(ierr);
>   }
>  }
>  ierr = MatAssemblyBegin(outMat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>  ierr = MatAssemblyEnd(outMat, MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);
>  return 0;
> }
>
>


-- 
What most experimenters take for granted before they begin their experiments
is infinitely more interesting than any results to which their experiments
lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100518/093bb220/attachment.htm>