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