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