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 If a Mat type does not implement the operation, each processor for which this is called 18 gets the values for its rows using MatGetRow(). 19 20 The vector must have the same parallel row layout as the matrix. 21 22 Contributed by: Denis Vanderstraeten 23 24 .seealso: MatGetRow(), MatGetDiagonal(), MatMult() 25 26 @*/ 27 PetscErrorCode MatGetColumnVector(Mat A,Vec yy,PetscInt col) 28 { 29 PetscScalar *y; 30 const PetscScalar *v; 31 PetscErrorCode ierr; 32 PetscInt i,j,nz,N,Rs,Re,rs,re; 33 const PetscInt *idx; 34 35 PetscFunctionBegin; 36 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 37 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 38 PetscValidLogicalCollectiveInt(A,col,3); 39 if (col < 0) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col); 40 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 41 if (col >= N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N); 42 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 43 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 44 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); 45 46 if (A->ops->getcolumnvector) { 47 ierr = (*A->ops->getcolumnvector)(A,yy,col);CHKERRQ(ierr); 48 } else { 49 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 50 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 51 /* TODO for general matrices */ 52 for (i=Rs; i<Re; i++) { 53 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 54 if (nz && idx[0] <= col) { 55 /* 56 Should use faster search here 57 */ 58 for (j=0; j<nz; j++) { 59 if (idx[j] >= col) { 60 if (idx[j] == col) y[i-rs] = v[j]; 61 break; 62 } 63 } 64 } 65 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 66 } 67 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /*@ 73 MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix. 74 75 Input Parameter: 76 + A - the matrix 77 - type - NORM_2, NORM_1 or NORM_INFINITY 78 79 Output Parameter: 80 . norms - an array as large as the TOTAL number of columns in the matrix 81 82 Level: intermediate 83 84 Notes: 85 Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values, 86 if each process wants only some of the values it should extract the ones it wants from the array. 87 88 .seealso: NormType, MatNorm() 89 90 @*/ 91 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[]) 92 { 93 PetscErrorCode ierr; 94 ReductionType reductiontype; 95 96 PetscFunctionBegin; 97 switch(type) { 98 case NORM_2: 99 reductiontype = REDUCTION_NORM_2; 100 break; 101 case NORM_1: 102 reductiontype = REDUCTION_NORM_1; 103 break; 104 case NORM_FROBENIUS: 105 reductiontype = REDUCTION_NORM_FROBENIUS; 106 break; 107 case NORM_INFINITY: 108 reductiontype = REDUCTION_NORM_INFINITY; 109 break; 110 case NORM_1_AND_2: 111 reductiontype = REDUCTION_NORM_1_AND_2; 112 break; 113 default: 114 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 115 } 116 ierr = MatGetColumnReductions(A,reductiontype,norms);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 /*@ 121 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 122 123 Input Parameter: 124 + A - the matrix 125 - type - NORM_2, NORM_1, NORM_INFINITY, SUM, MEAN 126 127 Output Parameter: 128 . reductions - an array as large as the TOTAL number of columns in the matrix 129 130 Level: intermediate 131 132 Notes: 133 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 134 if each process wants only some of the values it should extract the ones it wants from the array. 135 136 Developer Note: 137 MatGetColumnNorms() is now implemented using this routine. 138 139 .seealso: NormType, MatGetColumnNorms() 140 141 @*/ 142 PetscErrorCode MatGetColumnReductions(Mat A,ReductionType type,PetscReal reductions[]) 143 { 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 148 if (A->ops->getcolumnreductions) { 149 ierr = (*A->ops->getcolumnreductions)(A,type,reductions);CHKERRQ(ierr); 150 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type"); 151 PetscFunctionReturn(0); 152 } 153