Actual source code: getcolv.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>

  4: /*@
  5:    MatGetColumnVector - Gets the values from a given column of a matrix.

  7:    Not Collective

  9:    Input Parameters:
 10: +  A - the matrix
 11: .  yy - the vector
 12: -  col - the column requested (in global numbering)

 14:    Level: advanced

 16:    Notes:
 17:    Each processor for which this is called gets the values for its rows.

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

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

 24:    Contributed by: Denis Vanderstraeten

 26: .keywords: matrix, column, get

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

 30: @*/
 31: PetscErrorCode  MatGetColumnVector(Mat A,Vec yy,PetscInt col)
 32: {
 33:   PetscScalar       *y;
 34:   const PetscScalar *v;
 35:   PetscErrorCode    ierr;
 36:   PetscInt          i,j,nz,N,Rs,Re,rs,re;
 37:   const PetscInt    *idx;

 42:   if (col < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col);
 43:   MatGetSize(A,NULL,&N);
 44:   if (col >= N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N);
 45:   MatGetOwnershipRange(A,&Rs,&Re);
 46:   VecGetOwnershipRange(yy,&rs,&re);
 47:   if (Rs != rs || Re != re) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Matrix %D %D does not have same ownership range (size) as vector %D %D",Rs,Re,rs,re);

 49:   if (A->ops->getcolumnvector) {
 50:     (*A->ops->getcolumnvector)(A,yy,col);
 51:   } else {
 52:     VecSet(yy,0.0);
 53:     VecGetArray(yy,&y);

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




 78: /*@
 79:     MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.

 81:   Input Parameter:
 82: +  A - the matrix
 83: -  type - NORM_2, NORM_1 or NORM_INFINITY

 85:   Output Parameter:
 86: .  norms - an array as large as the TOTAL number of columns in the matrix

 88:    Level: intermediate

 90:    Notes:
 91:     Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values,
 92:     if each process wants only some of the values it should extract the ones it wants from the array.

 94: .seealso: NormType, MatNorm()

 96: @*/
 97: PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[])
 98: {

103:   if (A->ops->getcolumnnorms) {
104:     (*A->ops->getcolumnnorms)(A,type,norms);
105:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type");
106:   return(0);
107: }