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 PetscInt i,j,nz,N,Rs,Re,rs,re; 32 const PetscInt *idx; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 36 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 37 PetscValidLogicalCollectiveInt(A,col,3); 38 PetscCheck(col >= 0,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %" PetscInt_FMT,col); 39 PetscCall(MatGetSize(A,NULL,&N)); 40 PetscCheck(col < N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested column %" PetscInt_FMT " larger than number columns in matrix %" PetscInt_FMT,col,N); 41 PetscCall(MatGetOwnershipRange(A,&Rs,&Re)); 42 PetscCall(VecGetOwnershipRange(yy,&rs,&re)); 43 PetscCheck(Rs == rs && Re == re,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Matrix %" PetscInt_FMT " %" PetscInt_FMT " does not have same ownership range (size) as vector %" PetscInt_FMT " %" PetscInt_FMT,Rs,Re,rs,re); 44 45 if (A->ops->getcolumnvector) { 46 PetscCall((*A->ops->getcolumnvector)(A,yy,col)); 47 } else { 48 PetscCall(VecSet(yy,0.0)); 49 PetscCall(VecGetArray(yy,&y)); 50 /* TODO for general matrices */ 51 for (i=Rs; i<Re; i++) { 52 PetscCall(MatGetRow(A,i,&nz,&idx,&v)); 53 if (nz && idx[0] <= col) { 54 /* 55 Should use faster search here 56 */ 57 for (j=0; j<nz; j++) { 58 if (idx[j] >= col) { 59 if (idx[j] == col) y[i-rs] = v[j]; 60 break; 61 } 62 } 63 } 64 PetscCall(MatRestoreRow(A,i,&nz,&idx,&v)); 65 } 66 PetscCall(VecRestoreArray(yy,&y)); 67 } 68 PetscFunctionReturn(0); 69 } 70 71 /*@ 72 MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix. 73 74 Input Parameters: 75 + A - the matrix 76 - type - NORM_2, NORM_1 or NORM_INFINITY 77 78 Output Parameter: 79 . norms - an array as large as the TOTAL number of columns in the matrix 80 81 Level: intermediate 82 83 Notes: 84 Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values, 85 if each process wants only some of the values it should extract the ones it wants from the array. 86 87 .seealso: `NormType`, `MatNorm()` 88 89 @*/ 90 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[]) 91 { 92 /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions(). 93 * I've kept this as a function because it allows slightly more in the way of error checking, 94 * erroring out if MatGetColumnNorms() is not called with a valid NormType. */ 95 96 PetscFunctionBegin; 97 if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) { 98 PetscCall(MatGetColumnReductions(A,type,norms)); 99 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 100 PetscFunctionReturn(0); 101 } 102 103 /*@ 104 MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix. 105 106 Input Parameter: 107 . A - the matrix 108 109 Output Parameter: 110 . sums - an array as large as the TOTAL number of columns in the matrix 111 112 Level: intermediate 113 114 Notes: 115 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 116 if each process wants only some of the values it should extract the ones it wants from the array. 117 118 .seealso: `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 119 120 @*/ 121 PetscErrorCode MatGetColumnSumsRealPart(Mat A,PetscReal sums[]) 122 { 123 PetscFunctionBegin; 124 PetscCall(MatGetColumnReductions(A,REDUCTION_SUM_REALPART,sums)); 125 PetscFunctionReturn(0); 126 } 127 128 /*@ 129 MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix. 130 131 Input Parameter: 132 . A - the matrix 133 134 Output Parameter: 135 . sums - an array as large as the TOTAL number of columns in the matrix 136 137 Level: intermediate 138 139 Notes: 140 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 141 if each process wants only some of the values it should extract the ones it wants from the array. 142 143 .seealso: `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 144 145 @*/ 146 PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A,PetscReal sums[]) 147 { 148 PetscFunctionBegin; 149 PetscCall(MatGetColumnReductions(A,REDUCTION_SUM_IMAGINARYPART,sums)); 150 PetscFunctionReturn(0); 151 } 152 153 /*@ 154 MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix. 155 156 Input Parameter: 157 . A - the matrix 158 159 Output Parameter: 160 . sums - an array as large as the TOTAL number of columns in the matrix 161 162 Level: intermediate 163 164 Notes: 165 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 166 if each process wants only some of the values it should extract the ones it wants from the array. 167 168 .seealso: `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 169 170 @*/ 171 PetscErrorCode MatGetColumnSums(Mat A,PetscScalar sums[]) 172 { 173 #if defined(PETSC_USE_COMPLEX) 174 PetscInt i,n; 175 PetscReal *work; 176 #endif 177 178 PetscFunctionBegin; 179 180 #if !defined(PETSC_USE_COMPLEX) 181 PetscCall(MatGetColumnSumsRealPart(A,sums)); 182 #else 183 PetscCall(MatGetSize(A,NULL,&n)); 184 PetscCall(PetscArrayzero(sums,n)); 185 PetscCall(PetscCalloc1(n,&work)); 186 PetscCall(MatGetColumnSumsRealPart(A,work)); 187 for (i=0; i<n; i++) sums[i] = work[i]; 188 PetscCall(MatGetColumnSumsImaginaryPart(A,work)); 189 for (i=0; i<n; i++) sums[i] += work[i]*PETSC_i; 190 PetscCall(PetscFree(work)); 191 #endif 192 PetscFunctionReturn(0); 193 } 194 195 /*@ 196 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix. 197 198 Input Parameter: 199 . A - the matrix 200 201 Output Parameter: 202 . sums - an array as large as the TOTAL number of columns in the matrix 203 204 Level: intermediate 205 206 Notes: 207 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 208 if each process wants only some of the values it should extract the ones it wants from the array. 209 210 .seealso: `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 211 212 @*/ 213 PetscErrorCode MatGetColumnMeansRealPart(Mat A,PetscReal means[]) 214 { 215 PetscFunctionBegin; 216 PetscCall(MatGetColumnReductions(A,REDUCTION_MEAN_REALPART,means)); 217 PetscFunctionReturn(0); 218 } 219 220 /*@ 221 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix. 222 223 Input Parameter: 224 . A - the matrix 225 226 Output Parameter: 227 . sums - an array as large as the TOTAL number of columns in the matrix 228 229 Level: intermediate 230 231 Notes: 232 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 233 if each process wants only some of the values it should extract the ones it wants from the array. 234 235 .seealso: `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 236 237 @*/ 238 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A,PetscReal means[]) 239 { 240 PetscFunctionBegin; 241 PetscCall(MatGetColumnReductions(A,REDUCTION_MEAN_IMAGINARYPART,means)); 242 PetscFunctionReturn(0); 243 } 244 245 /*@ 246 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix. 247 248 Input Parameter: 249 . A - the matrix 250 251 Output Parameter: 252 . means - an array as large as the TOTAL number of columns in the matrix 253 254 Level: intermediate 255 256 Notes: 257 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 258 if each process wants only some of the values it should extract the ones it wants from the array. 259 260 .seealso: `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 261 262 @*/ 263 PetscErrorCode MatGetColumnMeans(Mat A,PetscScalar means[]) 264 { 265 #if defined(PETSC_USE_COMPLEX) 266 PetscInt i,n; 267 PetscReal *work; 268 #endif 269 270 PetscFunctionBegin; 271 272 #if !defined(PETSC_USE_COMPLEX) 273 PetscCall(MatGetColumnMeansRealPart(A,means)); 274 #else 275 PetscCall(MatGetSize(A,NULL,&n)); 276 PetscCall(PetscArrayzero(means,n)); 277 PetscCall(PetscCalloc1(n,&work)); 278 PetscCall(MatGetColumnMeansRealPart(A,work)); 279 for (i=0; i<n; i++) means[i] = work[i]; 280 PetscCall(MatGetColumnMeansImaginaryPart(A,work)); 281 for (i=0; i<n; i++) means[i] += work[i]*PETSC_i; 282 PetscCall(PetscFree(work)); 283 #endif 284 PetscFunctionReturn(0); 285 } 286 287 /*@ 288 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 289 290 Input Parameters: 291 + A - the matrix 292 - type - A constant defined in NormType or ReductionType: NORM_2, NORM_1, NORM_INFINITY, REDUCTION_SUM_REALPART, 293 REDUCTION_SUM_IMAGINARYPART, REDUCTION_MEAN_REALPART, REDUCTION_MEAN_IMAGINARYPART 294 295 Output Parameter: 296 . reductions - an array as large as the TOTAL number of columns in the matrix 297 298 Level: developer 299 300 Notes: 301 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 302 if each process wants only some of the values it should extract the ones it wants from the array. 303 304 Developer Note: 305 This routine is primarily intended as a back-end. 306 MatGetColumnNorms(), MatGetColumnSums(), and MatGetColumnMeans() are implemented using this routine. 307 308 .seealso: `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()` 309 310 @*/ 311 PetscErrorCode MatGetColumnReductions(Mat A,PetscInt type,PetscReal reductions[]) 312 { 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 315 if (A->ops->getcolumnreductions) { 316 PetscCall((*A->ops->getcolumnreductions)(A,type,reductions)); 317 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type"); 318 PetscFunctionReturn(0); 319 } 320