1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 3 /*@ 4 MatGetColumnVector - Gets the values from a given column of a matrix. 5 6 Not Collective 7 8 Input Parameters: 9 + A - the matrix 10 . yy - the vector 11 - col - the column requested (in global numbering) 12 13 Level: advanced 14 15 Notes: 16 If a `MatType` does not implement the operation, each processor for which this is called 17 gets the values for its rows using `MatGetRow()`. 18 19 The vector must have the same parallel row layout as the matrix. 20 21 Contributed by: Denis Vanderstraeten 22 23 .seealso: [](ch_matrices), `Mat`, `MatGetRow()`, `MatGetDiagonal()`, `MatMult()` 24 @*/ 25 PetscErrorCode MatGetColumnVector(Mat A, Vec yy, PetscInt col) 26 { 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(PETSC_SUCCESS); 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: [](ch_matrices), `Mat`, `NormType`, `MatNorm()` 85 @*/ 86 PetscErrorCode MatGetColumnNorms(Mat A, NormType type, PetscReal norms[]) 87 { 88 /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions(). 89 * I've kept this as a function because it allows slightly more in the way of error checking, 90 * erroring out if MatGetColumnNorms() is not called with a valid NormType. */ 91 92 PetscFunctionBegin; 93 if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) { 94 PetscCall(MatGetColumnReductions(A, type, norms)); 95 } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown NormType"); 96 PetscFunctionReturn(PETSC_SUCCESS); 97 } 98 99 /*@ 100 MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix. 101 102 Input Parameter: 103 . A - the matrix 104 105 Output Parameter: 106 . sums - an array as large as the TOTAL number of columns in the matrix 107 108 Level: intermediate 109 110 Note: 111 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 112 if each process wants only some of the values it should extract the ones it wants from the array. 113 114 .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 115 @*/ 116 PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[]) 117 { 118 PetscFunctionBegin; 119 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums)); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 /*@ 124 MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix. 125 126 Input Parameter: 127 . A - the matrix 128 129 Output Parameter: 130 . sums - an array as large as the TOTAL number of columns in the matrix 131 132 Level: intermediate 133 134 Note: 135 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 136 if each process wants only some of the values it should extract the ones it wants from the array. 137 138 .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 139 @*/ 140 PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[]) 141 { 142 PetscFunctionBegin; 143 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums)); 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 /*@ 148 MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix. 149 150 Input Parameter: 151 . A - the matrix 152 153 Output Parameter: 154 . sums - an array as large as the TOTAL number of columns in the matrix 155 156 Level: intermediate 157 158 Note: 159 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 160 if each process wants only some of the values it should extract the ones it wants from the array. 161 162 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 163 @*/ 164 PetscErrorCode MatGetColumnSums(Mat A, PetscScalar sums[]) 165 { 166 #if defined(PETSC_USE_COMPLEX) 167 PetscInt i, n; 168 PetscReal *work; 169 #endif 170 171 PetscFunctionBegin; 172 173 #if !defined(PETSC_USE_COMPLEX) 174 PetscCall(MatGetColumnSumsRealPart(A, sums)); 175 #else 176 PetscCall(MatGetSize(A, NULL, &n)); 177 PetscCall(PetscArrayzero(sums, n)); 178 PetscCall(PetscCalloc1(n, &work)); 179 PetscCall(MatGetColumnSumsRealPart(A, work)); 180 for (i = 0; i < n; i++) sums[i] = work[i]; 181 PetscCall(MatGetColumnSumsImaginaryPart(A, work)); 182 for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i; 183 PetscCall(PetscFree(work)); 184 #endif 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 /*@ 189 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix. 190 191 Input Parameter: 192 . A - the matrix 193 194 Output Parameter: 195 . means - an array as large as the TOTAL number of columns in the matrix 196 197 Level: intermediate 198 199 Note: 200 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 201 if each process wants only some of the values it should extract the ones it wants from the array. 202 203 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 204 @*/ 205 PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[]) 206 { 207 PetscFunctionBegin; 208 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means)); 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 /*@ 213 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix. 214 215 Input Parameter: 216 . A - the matrix 217 218 Output Parameter: 219 . means - an array as large as the TOTAL number of columns in the matrix 220 221 Level: intermediate 222 223 Note: 224 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 225 if each process wants only some of the values it should extract the ones it wants from the array. 226 227 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 228 @*/ 229 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[]) 230 { 231 PetscFunctionBegin; 232 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means)); 233 PetscFunctionReturn(PETSC_SUCCESS); 234 } 235 236 /*@ 237 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix. 238 239 Input Parameter: 240 . A - the matrix 241 242 Output Parameter: 243 . means - an array as large as the TOTAL number of columns in the matrix 244 245 Level: intermediate 246 247 Note: 248 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 249 if each process wants only some of the values it should extract the ones it wants from the array. 250 251 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 252 @*/ 253 PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[]) 254 { 255 #if defined(PETSC_USE_COMPLEX) 256 PetscInt i, n; 257 PetscReal *work; 258 #endif 259 260 PetscFunctionBegin; 261 262 #if !defined(PETSC_USE_COMPLEX) 263 PetscCall(MatGetColumnMeansRealPart(A, means)); 264 #else 265 PetscCall(MatGetSize(A, NULL, &n)); 266 PetscCall(PetscArrayzero(means, n)); 267 PetscCall(PetscCalloc1(n, &work)); 268 PetscCall(MatGetColumnMeansRealPart(A, work)); 269 for (i = 0; i < n; i++) means[i] = work[i]; 270 PetscCall(MatGetColumnMeansImaginaryPart(A, work)); 271 for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i; 272 PetscCall(PetscFree(work)); 273 #endif 274 PetscFunctionReturn(PETSC_SUCCESS); 275 } 276 277 /*@ 278 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 279 280 Input Parameters: 281 + A - the matrix 282 - type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`, 283 `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART` 284 285 Output Parameter: 286 . reductions - an array as large as the TOTAL number of columns in the matrix 287 288 Level: developer 289 290 Note: 291 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 292 if each process wants only some of the values it should extract the ones it wants from the array. 293 294 Developer Notes: 295 This routine is primarily intended as a back-end. 296 `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine. 297 298 .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()` 299 @*/ 300 PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[]) 301 { 302 PetscFunctionBegin; 303 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 304 PetscUseTypeMethod(A, getcolumnreductions, type, reductions); 305 PetscFunctionReturn(PETSC_SUCCESS); 306 } 307