1 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/ 2 3 typedef struct { 4 Mat A; 5 Mat D; /* local submatrix for diagonal part */ 6 Vec w; 7 } Mat_NormalHermitian; 8 9 static PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 10 { 11 Mat_NormalHermitian *a; 12 Mat B, *suba; 13 IS *row; 14 PetscScalar shift, scale; 15 PetscInt M; 16 17 PetscFunctionBegin; 18 PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented"); 19 PetscCall(MatShellGetScalingShifts(mat, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 20 PetscCall(MatShellGetContext(mat, &a)); 21 B = a->A; 22 if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat)); 23 PetscCall(MatGetSize(B, &M, NULL)); 24 PetscCall(PetscMalloc1(n, &row)); 25 PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0])); 26 PetscCall(ISSetIdentity(row[0])); 27 for (M = 1; M < n; ++M) row[M] = row[0]; 28 PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba)); 29 for (M = 0; M < n; ++M) { 30 PetscCall(MatCreateNormalHermitian(suba[M], *submat + M)); 31 PetscCall(MatShift((*submat)[M], shift)); 32 PetscCall(MatScale((*submat)[M], scale)); 33 } 34 PetscCall(ISDestroy(&row[0])); 35 PetscCall(PetscFree(row)); 36 PetscCall(MatDestroySubMatrices(n, &suba)); 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B) 41 { 42 Mat_NormalHermitian *a; 43 Mat C, Aa; 44 IS row; 45 PetscScalar shift, scale; 46 47 PetscFunctionBegin; 48 PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same"); 49 PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 50 PetscCall(MatShellGetContext(A, &a)); 51 Aa = a->A; 52 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row)); 53 PetscCall(ISSetIdentity(row)); 54 PetscCall(MatPermute(Aa, row, colp, &C)); 55 PetscCall(ISDestroy(&row)); 56 PetscCall(MatCreateNormalHermitian(C, B)); 57 PetscCall(MatDestroy(&C)); 58 PetscCall(MatShift(*B, shift)); 59 PetscCall(MatScale(*B, scale)); 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B) 64 { 65 Mat_NormalHermitian *a; 66 Mat C; 67 68 PetscFunctionBegin; 69 PetscCall(MatShellGetContext(A, &a)); 70 PetscCall(MatDuplicate(a->A, op, &C)); 71 PetscCall(MatCreateNormalHermitian(C, B)); 72 PetscCall(MatDestroy(&C)); 73 if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN)); 74 PetscFunctionReturn(PETSC_SUCCESS); 75 } 76 77 static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str) 78 { 79 Mat_NormalHermitian *a, *b; 80 81 PetscFunctionBegin; 82 PetscCall(MatShellGetContext(A, &a)); 83 PetscCall(MatShellGetContext(B, &b)); 84 PetscCall(MatCopy(a->A, b->A, str)); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y) 89 { 90 Mat_NormalHermitian *Na; 91 92 PetscFunctionBegin; 93 PetscCall(MatShellGetContext(N, &Na)); 94 PetscCall(MatMult(Na->A, x, Na->w)); 95 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y)); 96 PetscFunctionReturn(PETSC_SUCCESS); 97 } 98 99 static PetscErrorCode MatDestroy_NormalHermitian(Mat N) 100 { 101 Mat_NormalHermitian *Na; 102 103 PetscFunctionBegin; 104 PetscCall(MatShellGetContext(N, &Na)); 105 PetscCall(MatDestroy(&Na->A)); 106 PetscCall(MatDestroy(&Na->D)); 107 PetscCall(VecDestroy(&Na->w)); 108 PetscCall(PetscFree(Na)); 109 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalHermitianGetMat_C", NULL)); 110 #if !defined(PETSC_USE_COMPLEX) 111 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL)); 112 #endif 113 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL)); 114 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL)); 115 #if defined(PETSC_HAVE_HYPRE) 116 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL)); 117 #endif 118 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL)); 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 /* 123 Slow, nonscalable version 124 */ 125 static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v) 126 { 127 Mat_NormalHermitian *Na; 128 Mat A; 129 PetscInt i, j, rstart, rend, nnz; 130 const PetscInt *cols; 131 PetscScalar *diag, *work, *values; 132 const PetscScalar *mvalues; 133 PetscMPIInt iN; 134 135 PetscFunctionBegin; 136 PetscCall(MatShellGetContext(N, &Na)); 137 A = Na->A; 138 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work)); 139 PetscCall(PetscArrayzero(work, A->cmap->N)); 140 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 141 for (i = rstart; i < rend; i++) { 142 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues)); 143 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]); 144 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues)); 145 } 146 PetscCall(PetscMPIIntCast(A->cmap->N, &iN)); 147 PetscCallMPI(MPIU_Allreduce(work, diag, iN, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N))); 148 rstart = N->cmap->rstart; 149 rend = N->cmap->rend; 150 PetscCall(VecGetArray(v, &values)); 151 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart)); 152 PetscCall(VecRestoreArray(v, &values)); 153 PetscCall(PetscFree2(diag, work)); 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 157 static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D) 158 { 159 Mat_NormalHermitian *Na; 160 Mat M, A; 161 162 PetscFunctionBegin; 163 PetscCall(MatShellGetContext(N, &Na)); 164 A = Na->A; 165 PetscCall(MatGetDiagonalBlock(A, &M)); 166 PetscCall(MatCreateNormalHermitian(M, &Na->D)); 167 *D = Na->D; 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M) 172 { 173 Mat_NormalHermitian *Aa; 174 175 PetscFunctionBegin; 176 PetscCall(MatShellGetContext(A, &Aa)); 177 *M = Aa->A; 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 /*@ 182 MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN` 183 184 Logically Collective 185 186 Input Parameter: 187 . A - the `MATNORMALHERMITIAN` matrix 188 189 Output Parameter: 190 . M - the matrix object stored inside `A` 191 192 Level: intermediate 193 194 .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()` 195 @*/ 196 PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M) 197 { 198 PetscFunctionBegin; 199 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 200 PetscValidType(A, 1); 201 PetscAssertPointer(M, 2); 202 PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M)); 203 PetscFunctionReturn(PETSC_SUCCESS); 204 } 205 206 static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 207 { 208 Mat_NormalHermitian *Aa; 209 Mat B, conjugate; 210 PetscInt m, n, M, N; 211 212 PetscFunctionBegin; 213 PetscCall(MatShellGetContext(A, &Aa)); 214 PetscCall(MatGetSize(A, &M, &N)); 215 PetscCall(MatGetLocalSize(A, &m, &n)); 216 if (reuse == MAT_REUSE_MATRIX) { 217 B = *newmat; 218 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B)); 219 } else { 220 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B)); 221 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 222 PetscCall(MatProductSetFromOptions(B)); 223 PetscCall(MatProductSymbolic(B)); 224 PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE)); 225 } 226 if (PetscDefined(USE_COMPLEX)) { 227 PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate)); 228 PetscCall(MatConjugate(conjugate)); 229 PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B)); 230 } 231 PetscCall(MatProductNumeric(B)); 232 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate)); 233 if (reuse == MAT_INPLACE_MATRIX) { 234 PetscCall(MatHeaderReplace(A, &B)); 235 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 236 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat)); 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 #if defined(PETSC_HAVE_HYPRE) 241 static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B) 242 { 243 PetscFunctionBegin; 244 if (reuse == MAT_INITIAL_MATRIX) { 245 PetscCall(MatConvert(A, MATAIJ, reuse, B)); 246 PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B)); 247 } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */ 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 #endif 251 252 /*MC 253 MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A 254 255 Level: intermediate 256 257 Developer Notes: 258 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code 259 260 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage 261 262 .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()` 263 M*/ 264 265 /*@ 266 MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A. 267 268 Collective 269 270 Input Parameter: 271 . A - the (possibly rectangular complex) matrix 272 273 Output Parameter: 274 . N - the matrix that represents (A*)'*A 275 276 Level: intermediate 277 278 Note: 279 The product (A*)'*A is NOT actually formed! Rather the new matrix 280 object performs the matrix-vector product, `MatMult()`, by first multiplying by 281 A and then (A*)' 282 283 .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()` 284 @*/ 285 PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N) 286 { 287 Mat_NormalHermitian *Na; 288 VecType vtype; 289 290 PetscFunctionBegin; 291 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 292 PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap)); 293 PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap)); 294 PetscCall(MatSetType(*N, MATSHELL)); 295 PetscCall(PetscNew(&Na)); 296 PetscCall(MatShellSetContext(*N, Na)); 297 PetscCall(PetscObjectReference((PetscObject)A)); 298 Na->A = A; 299 PetscCall(MatCreateVecs(A, NULL, &Na->w)); 300 301 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 302 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_NormalHermitian)); 303 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_NormalHermitian)); 304 PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian)); 305 #if !defined(PETSC_USE_COMPLEX) 306 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian)); 307 #endif 308 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_NormalHermitian)); 309 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NormalHermitian)); 310 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (void (*)(void))MatGetDiagonalBlock_NormalHermitian)); 311 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_NormalHermitian)); 312 (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian; 313 (*N)->ops->permute = MatPermute_NormalHermitian; 314 315 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian)); 316 #if !defined(PETSC_USE_COMPLEX) 317 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalHermitianGetMat_NormalHermitian)); 318 #endif 319 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ)); 320 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ)); 321 #if defined(PETSC_HAVE_HYPRE) 322 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE)); 323 #endif 324 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable)); 325 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable)); 326 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable)); 327 PetscCall(MatSetOption(*N, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE)); 328 PetscCall(MatGetVecType(A, &vtype)); 329 PetscCall(MatSetVecType(*N, vtype)); 330 #if defined(PETSC_HAVE_DEVICE) 331 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 332 #endif 333 PetscCall(MatSetUp(*N)); 334 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN)); 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337