1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 /*@ 5 MatGetColumnVector - Gets the values from a given column of a matrix. 6 7 Not Collective 8 9 Input Parameters: 10 + A - the matrix 11 . yy - the vector 12 - col - the column requested (in global numbering) 13 14 Level: advanced 15 16 Notes: 17 Each processor for which this is called gets the values for its rows. 18 19 Since PETSc matrices are usually stored in compressed row format, this routine 20 will generally be slow. 21 22 The vector must have the same parallel row layout as the matrix. 23 24 Contributed by: Denis Vanderstraeten 25 26 .keywords: matrix, column, get 27 28 .seealso: MatGetRow(), MatGetDiagonal() 29 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; 38 39 PetscFunctionBegin; 40 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 41 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 42 if (col < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col); 43 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 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 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 46 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 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); 48 49 if (A->ops->getcolumnvector) { 50 ierr = (*A->ops->getcolumnvector)(A,yy,col);CHKERRQ(ierr); 51 } else { 52 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 53 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 54 55 for (i=Rs; i<Re; i++) { 56 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 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 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 69 } 70 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 71 } 72 PetscFunctionReturn(0); 73 } 74 75 76 77 78 /*@ 79 MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix. 80 81 Input Parameter: 82 + A - the matrix 83 - type - NORM_2, NORM_1 or NORM_INFINITY 84 85 Output Parameter: 86 . norms - an array as large as the TOTAL number of columns in the matrix 87 88 Level: intermediate 89 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. 93 94 .seealso: NormType, MatNorm() 95 96 @*/ 97 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[]) 98 { 99 PetscErrorCode ierr; 100 101 PetscFunctionBegin; 102 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 103 if (A->ops->getcolumnnorms) { 104 ierr = (*A->ops->getcolumnnorms)(A,type,norms);CHKERRQ(ierr); 105 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type"); 106 PetscFunctionReturn(0); 107 } 108 109