Actual source code: getcolv.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/matimpl.h

      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: PetscErrorCode PETSCMAT_DLLEXPORT MatGetColumnVector(Mat A,Vec yy,PetscInt col)
     35: {
     36:   PetscScalar        *y,zero = 0.0;
     37:   const PetscScalar  *v;
     38:   PetscErrorCode     ierr;
     39:   PetscInt           i,j,nz,N,Rs,Re,rs,re;
     40:   const PetscInt     *idx;
     41:   MPI_Comm           comm;
     42:

     47:   if (col < 0)  SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col);
     48:   MatGetSize(A,PETSC_NULL,&N);
     49:   if (col >= N)  SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N);

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

     53:   PetscObjectGetComm((PetscObject)yy,&comm);
     54:   VecGetOwnershipRange(yy,&rs,&re);
     55:   if (Rs != rs || Re != re) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Matrix %D %D does not have same ownership range (size) as vector %D %D",Rs,Re,rs,re);

     57:   VecSet(yy,zero);
     58:   VecGetArray(yy,&y);

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

     76:   VecRestoreArray(yy,&y);
     77:   return(0);
     78: }