1: /*$Id: getcolv.c,v 1.23 2001/08/08 15:47:19 bsmith Exp $*/

 3:  #include src/mat/matimpl.h

  5: #undef __FUNCT__  
  7: /*@
  8:    MatGetColumnVector - Gets the values from a given column of a matrix.

 10:    Not Collective

 12:    Input Parameters:
 13: +  A - the matrix
 14: .  yy - the vector
 15: -  c - the column requested (in global numbering)

 17:    Level: advanced

 19:    Notes:
 20:    Each processor for which this is called gets the values for its rows.

 22:    Since PETSc matrices are usually stored in compressed row format, this routine
 23:    will generally be slow.

 25:    The vector must have the same parallel row layout as the matrix.

 27:    Contributed by: Denis Vanderstraeten

 29: .keywords: matrix, column, get 

 31: .seealso: MatGetRow(), MatGetDiagonal()

 33: @*/
 34: int MatGetColumnVector(Mat A,Vec yy,int col)
 35: {
 36:   PetscScalar   *y,*v,zero = 0.0;
 37:   int      ierr,i,j,nz,*idx,N,Rs,Re,rs,re;
 38:   MPI_Comm comm;
 39: 

 44:   if (col < 0)  SETERRQ1(1,"Requested negative column: %d",col);
 45:   MatGetSize(A,PETSC_NULL,&N);
 46:   if (col >= N)  SETERRQ2(1,"Requested column %d larger than number columns in matrix %d",col,N);

 48:   MatGetOwnershipRange(A,&Rs,&Re);

 50:   PetscObjectGetComm((PetscObject)yy,&comm);
 51:   VecGetOwnershipRange(yy,&rs,&re);
 52:   if (Rs != rs || Re != re) SETERRQ4(1,"Matrix %d %d does not have same ownership range (size) as vector %d %d",Rs,Re,rs,re);

 54:   VecSet(&zero,yy);
 55:   VecGetArray(yy,&y);

 57:   for (i=Rs; i<Re; i++) {
 58:     MatGetRow(A,i,&nz,&idx,&v);
 59:     if (nz && idx[0] <= col) {
 60:       /*
 61:           Should use faster search here 
 62:       */
 63:       for (j=0; j<nz; j++) {
 64:         if (idx[j] >= col) {
 65:           if (idx[j] == col) y[i-rs] = v[j];
 66:           break;
 67:         }
 68:       }
 69:     }
 70:     MatRestoreRow(A,i,&nz,&idx,&v);
 71:   }

 73:   VecRestoreArray(yy,&y);
 74:   return(0);
 75: }