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 PetscCheck(type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown NormType"); 94 PetscCall(MatGetColumnReductions(A, type, norms)); 95 PetscFunctionReturn(PETSC_SUCCESS); 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: [](ch_matrices), `Mat`, `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 114 @*/ 115 PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[]) 116 { 117 PetscFunctionBegin; 118 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums)); 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 /*@ 123 MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix. 124 125 Input Parameter: 126 . A - the matrix 127 128 Output Parameter: 129 . sums - an array as large as the TOTAL number of columns in the matrix 130 131 Level: intermediate 132 133 Note: 134 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 135 if each process wants only some of the values it should extract the ones it wants from the array. 136 137 .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 138 @*/ 139 PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[]) 140 { 141 PetscFunctionBegin; 142 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /*@ 147 MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix. 148 149 Input Parameter: 150 . A - the matrix 151 152 Output Parameter: 153 . sums - an array as large as the TOTAL number of columns in the matrix 154 155 Level: intermediate 156 157 Note: 158 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 159 if each process wants only some of the values it should extract the ones it wants from the array. 160 161 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 162 @*/ 163 PetscErrorCode MatGetColumnSums(Mat A, PetscScalar sums[]) 164 { 165 #if defined(PETSC_USE_COMPLEX) 166 PetscInt i, n; 167 PetscReal *work; 168 #endif 169 170 PetscFunctionBegin; 171 #if !defined(PETSC_USE_COMPLEX) 172 PetscCall(MatGetColumnSumsRealPart(A, sums)); 173 #else 174 PetscCall(MatGetSize(A, NULL, &n)); 175 PetscCall(PetscArrayzero(sums, n)); 176 PetscCall(PetscCalloc1(n, &work)); 177 PetscCall(MatGetColumnSumsRealPart(A, work)); 178 for (i = 0; i < n; i++) sums[i] = work[i]; 179 PetscCall(MatGetColumnSumsImaginaryPart(A, work)); 180 for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i; 181 PetscCall(PetscFree(work)); 182 #endif 183 PetscFunctionReturn(PETSC_SUCCESS); 184 } 185 186 /*@ 187 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix. 188 189 Input Parameter: 190 . A - the matrix 191 192 Output Parameter: 193 . means - an array as large as the TOTAL number of columns in the matrix 194 195 Level: intermediate 196 197 Note: 198 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 199 if each process wants only some of the values it should extract the ones it wants from the array. 200 201 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 202 @*/ 203 PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[]) 204 { 205 PetscFunctionBegin; 206 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means)); 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 /*@ 211 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix. 212 213 Input Parameter: 214 . A - the matrix 215 216 Output Parameter: 217 . means - an array as large as the TOTAL number of columns in the matrix 218 219 Level: intermediate 220 221 Note: 222 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 223 if each process wants only some of the values it should extract the ones it wants from the array. 224 225 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 226 @*/ 227 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[]) 228 { 229 PetscFunctionBegin; 230 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means)); 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 234 /*@ 235 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix. 236 237 Input Parameter: 238 . A - the matrix 239 240 Output Parameter: 241 . means - an array as large as the TOTAL number of columns in the matrix 242 243 Level: intermediate 244 245 Note: 246 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 247 if each process wants only some of the values it should extract the ones it wants from the array. 248 249 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()` 250 @*/ 251 PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[]) 252 { 253 #if defined(PETSC_USE_COMPLEX) 254 PetscInt i, n; 255 PetscReal *work; 256 #endif 257 258 PetscFunctionBegin; 259 #if !defined(PETSC_USE_COMPLEX) 260 PetscCall(MatGetColumnMeansRealPart(A, means)); 261 #else 262 PetscCall(MatGetSize(A, NULL, &n)); 263 PetscCall(PetscArrayzero(means, n)); 264 PetscCall(PetscCalloc1(n, &work)); 265 PetscCall(MatGetColumnMeansRealPart(A, work)); 266 for (i = 0; i < n; i++) means[i] = work[i]; 267 PetscCall(MatGetColumnMeansImaginaryPart(A, work)); 268 for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i; 269 PetscCall(PetscFree(work)); 270 #endif 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 /*@ 275 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 276 277 Input Parameters: 278 + A - the matrix 279 - type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`, 280 `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART` 281 282 Output Parameter: 283 . reductions - an array as large as the TOTAL number of columns in the matrix 284 285 Level: developer 286 287 Note: 288 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 289 if each process wants only some of the values it should extract the ones it wants from the array. 290 291 Developer Notes: 292 This routine is primarily intended as a back-end. 293 `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine. 294 295 .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()` 296 @*/ 297 PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[]) 298 { 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 301 PetscUseTypeMethod(A, getcolumnreductions, type, reductions); 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304