2: #include <petsc-private/matimpl.h>  /*I   "petscmat.h"  I*/

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

  9:    Not Collective

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

 16:    Level: advanced

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

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

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

 26:    Contributed by: Denis Vanderstraeten

 28: .keywords: matrix, column, get 

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

 32: @*/
 33: PetscErrorCode  MatGetColumnVector(Mat A,Vec yy,PetscInt col)
 34: {
 35:   PetscScalar        *y;
 36:   const PetscScalar  *v;
 37:   PetscErrorCode     ierr;
 38:   PetscInt           i,j,nz,N,Rs,Re,rs,re;
 39:   const PetscInt     *idx;
 40: 
 44:   if (col < 0)  SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col);
 45:   MatGetSize(A,PETSC_NULL,&N);
 46:   if (col >= N)  SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N);
 47:   MatGetOwnershipRange(A,&Rs,&Re);
 48:   VecGetOwnershipRange(yy,&rs,&re);
 49:   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);

 51:   if (A->ops->getcolumnvector) {
 52:     (*A->ops->getcolumnvector)(A,yy,col);
 53:   } else {
 54:     VecSet(yy,0.0);
 55:     VecGetArray(yy,&y);
 56: 
 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:     }
 72:     VecRestoreArray(yy,&y);
 73:   }
 74:   return(0);
 75: }




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

 85:   Input Parameter:
 86: +  A - the matrix
 87: -  type - NORM_2, NORM_1 or NORM_INFINITY

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

 92:    Level: intermediate

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

 97: .seealso: MatGetColumns()

 99: @*/
100: PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal *norms)
101: {

106:   if (A->ops->getcolumnnorms) {
107:     (*A->ops->getcolumnnorms)(A,type,norms);
108:   } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not coded for this matrix type");
109:   return(0);
110: }
111: