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
MatCreateSubMatrices_NormalHermitian(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat * submat[])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
MatPermute_NormalHermitian(Mat A,IS rowp,IS colp,Mat * B)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
MatDuplicate_NormalHermitian(Mat A,MatDuplicateOption op,Mat * B)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
MatCopy_NormalHermitian(Mat A,Mat B,MatStructure str)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
MatMult_NormalHermitian(Mat N,Vec x,Vec y)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
MatDestroy_NormalHermitian(Mat N)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 */
MatGetDiagonal_NormalHermitian(Mat N,Vec v)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 *work, *values;
132 const PetscScalar *mvalues;
133
134 PetscFunctionBegin;
135 PetscCall(MatShellGetContext(N, &Na));
136 A = Na->A;
137 PetscCall(PetscMalloc1(A->cmap->N, &work));
138 PetscCall(PetscArrayzero(work, A->cmap->N));
139 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
140 for (i = rstart; i < rend; i++) {
141 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
142 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
143 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
144 }
145 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, work, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
146 rstart = N->cmap->rstart;
147 rend = N->cmap->rend;
148 PetscCall(VecGetArray(v, &values));
149 PetscCall(PetscArraycpy(values, work + rstart, rend - rstart));
150 PetscCall(VecRestoreArray(v, &values));
151 PetscCall(PetscFree(work));
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
MatGetDiagonalBlock_NormalHermitian(Mat N,Mat * D)155 static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
156 {
157 Mat_NormalHermitian *Na;
158 Mat M, A;
159
160 PetscFunctionBegin;
161 PetscCall(MatShellGetContext(N, &Na));
162 A = Na->A;
163 PetscCall(MatGetDiagonalBlock(A, &M));
164 PetscCall(MatCreateNormalHermitian(M, &Na->D));
165 *D = Na->D;
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
MatNormalHermitianGetMat_NormalHermitian(Mat A,Mat * M)169 static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
170 {
171 Mat_NormalHermitian *Aa;
172
173 PetscFunctionBegin;
174 PetscCall(MatShellGetContext(A, &Aa));
175 *M = Aa->A;
176 PetscFunctionReturn(PETSC_SUCCESS);
177 }
178
179 /*@
180 MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`
181
182 Logically Collective
183
184 Input Parameter:
185 . A - the `MATNORMALHERMITIAN` matrix
186
187 Output Parameter:
188 . M - the matrix object stored inside `A`
189
190 Level: intermediate
191
192 .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
193 @*/
MatNormalHermitianGetMat(Mat A,Mat * M)194 PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
195 {
196 PetscFunctionBegin;
197 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
198 PetscValidType(A, 1);
199 PetscAssertPointer(M, 2);
200 PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
MatConvert_NormalHermitian_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)204 static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
205 {
206 Mat_NormalHermitian *Aa;
207 Mat B, conjugate;
208 Vec left, right, dshift;
209 PetscScalar scale, shift;
210 PetscInt m, n, M, N;
211
212 PetscFunctionBegin;
213 PetscCall(MatShellGetContext(A, &Aa));
214 PetscCall(MatShellGetScalingShifts(A, &shift, &scale, &dshift, &left, &right, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
215 PetscCall(MatGetSize(A, &M, &N));
216 PetscCall(MatGetLocalSize(A, &m, &n));
217 if (reuse == MAT_REUSE_MATRIX) {
218 B = *newmat;
219 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
220 } else {
221 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
222 PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
223 PetscCall(MatProductSetFromOptions(B));
224 PetscCall(MatProductSymbolic(B));
225 PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
226 }
227 if (PetscDefined(USE_COMPLEX)) {
228 PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
229 PetscCall(MatConjugate(conjugate));
230 PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
231 }
232 PetscCall(MatProductNumeric(B));
233 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
234 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
235 else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
236 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
237 PetscCall(MatDiagonalScale(*newmat, left, right));
238 PetscCall(MatScale(*newmat, scale));
239 PetscCall(MatShift(*newmat, shift));
240 if (dshift) PetscCall(MatDiagonalSet(*newmat, dshift, ADD_VALUES));
241 PetscFunctionReturn(PETSC_SUCCESS);
242 }
243
244 #if defined(PETSC_HAVE_HYPRE)
MatConvert_NormalHermitian_HYPRE(Mat A,MatType type,MatReuse reuse,Mat * B)245 static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
246 {
247 PetscFunctionBegin;
248 if (reuse == MAT_INITIAL_MATRIX) {
249 PetscCall(MatConvert(A, MATAIJ, reuse, B));
250 PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
251 } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
252 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254 #endif
255
256 /*MC
257 MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A
258
259 Level: intermediate
260
261 Developer Notes:
262 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
263
264 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
265
266 .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
267 M*/
268
269 /*@
270 MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like $A^* A$.
271
272 Collective
273
274 Input Parameter:
275 . A - the (possibly rectangular complex) matrix
276
277 Output Parameter:
278 . N - the matrix that represents $ A^* A$
279
280 Level: intermediate
281
282 Note:
283 The product $ A^* A$ is NOT actually formed! Rather the new matrix
284 object performs the matrix-vector product, `MatMult()`, by first multiplying by
285 $A$ and then $A^*$
286
287 If `MatGetFactor()` is called on this matrix with `MAT_FACTOR_QR` then the inner matrix `A` is used for the factorization
288
289 .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
290 @*/
MatCreateNormalHermitian(Mat A,Mat * N)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(MatSetBlockSize(*N, A->cmap->bs));
308 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_NormalHermitian));
309 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_NormalHermitian));
310 PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMult_NormalHermitian));
311 #if !defined(PETSC_USE_COMPLEX)
312 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMult_NormalHermitian));
313 #endif
314 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_NormalHermitian));
315 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_NormalHermitian));
316 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (PetscErrorCodeFn *)MatGetDiagonalBlock_NormalHermitian));
317 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_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, 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