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