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