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 #if !defined(PETSC_USE_COMPLEX) 173 PetscCall(MatGetColumnSumsRealPart(A, sums)); 174 #else 175 PetscCall(MatGetSize(A, NULL, &n)); 176 PetscCall(PetscArrayzero(sums, n)); 177 PetscCall(PetscCalloc1(n, &work)); 178 PetscCall(MatGetColumnSumsRealPart(A, work)); 179 for (i = 0; i < n; i++) sums[i] = work[i]; 180 PetscCall(MatGetColumnSumsImaginaryPart(A, work)); 181 for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i; 182 PetscCall(PetscFree(work)); 183 #endif 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 /*@ 188 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix. 189 190 Input Parameter: 191 . A - the matrix 192 193 Output Parameter: 194 . means - an array as large as the TOTAL number of columns in the matrix 195 196 Level: intermediate 197 198 Note: 199 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 200 if each process wants only some of the values it should extract the ones it wants from the array. 201 202 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 203 @*/ 204 PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[]) 205 { 206 PetscFunctionBegin; 207 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means)); 208 PetscFunctionReturn(PETSC_SUCCESS); 209 } 210 211 /*@ 212 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix. 213 214 Input Parameter: 215 . A - the matrix 216 217 Output Parameter: 218 . means - an array as large as the TOTAL number of columns in the matrix 219 220 Level: intermediate 221 222 Note: 223 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 224 if each process wants only some of the values it should extract the ones it wants from the array. 225 226 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 227 @*/ 228 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[]) 229 { 230 PetscFunctionBegin; 231 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means)); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 /*@ 236 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix. 237 238 Input Parameter: 239 . A - the matrix 240 241 Output Parameter: 242 . means - an array as large as the TOTAL number of columns in the matrix 243 244 Level: intermediate 245 246 Note: 247 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 248 if each process wants only some of the values it should extract the ones it wants from the array. 249 250 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 251 @*/ 252 PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[]) 253 { 254 #if defined(PETSC_USE_COMPLEX) 255 PetscInt i, n; 256 PetscReal *work; 257 #endif 258 259 PetscFunctionBegin; 260 #if !defined(PETSC_USE_COMPLEX) 261 PetscCall(MatGetColumnMeansRealPart(A, means)); 262 #else 263 PetscCall(MatGetSize(A, NULL, &n)); 264 PetscCall(PetscArrayzero(means, n)); 265 PetscCall(PetscCalloc1(n, &work)); 266 PetscCall(MatGetColumnMeansRealPart(A, work)); 267 for (i = 0; i < n; i++) means[i] = work[i]; 268 PetscCall(MatGetColumnMeansImaginaryPart(A, work)); 269 for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i; 270 PetscCall(PetscFree(work)); 271 #endif 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 /*@ 276 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 277 278 Input Parameters: 279 + A - the matrix 280 - type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`, 281 `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART` 282 283 Output Parameter: 284 . reductions - an array as large as the TOTAL number of columns in the matrix 285 286 Level: developer 287 288 Note: 289 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 290 if each process wants only some of the values it should extract the ones it wants from the array. 291 292 Developer Notes: 293 This routine is primarily intended as a back-end. 294 `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine. 295 296 .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()` 297 @*/ 298 PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[]) 299 { 300 PetscFunctionBegin; 301 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 302 PetscUseTypeMethod(A, getcolumnreductions, type, reductions); 303 PetscFunctionReturn(PETSC_SUCCESS); 304 } 305