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

[Umut Tabak]

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;
}