1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 Mat D; /* local submatrix for diagonal part */ 7 Vec w, left, right, leftwork, rightwork; 8 PetscScalar scale; 9 } Mat_Normal; 10 11 PetscErrorCode MatScale_NormalHermitian(Mat inA, PetscScalar scale) { 12 Mat_Normal *a = (Mat_Normal *)inA->data; 13 14 PetscFunctionBegin; 15 a->scale *= scale; 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode MatDiagonalScale_NormalHermitian(Mat inA, Vec left, Vec right) { 20 Mat_Normal *a = (Mat_Normal *)inA->data; 21 22 PetscFunctionBegin; 23 if (left) { 24 if (!a->left) { 25 PetscCall(VecDuplicate(left, &a->left)); 26 PetscCall(VecCopy(left, a->left)); 27 } else { 28 PetscCall(VecPointwiseMult(a->left, left, a->left)); 29 } 30 } 31 if (right) { 32 if (!a->right) { 33 PetscCall(VecDuplicate(right, &a->right)); 34 PetscCall(VecCopy(right, a->right)); 35 } else { 36 PetscCall(VecPointwiseMult(a->right, right, a->right)); 37 } 38 } 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) { 43 Mat_Normal *a = (Mat_Normal *)mat->data; 44 Mat B = a->A, *suba; 45 IS *row; 46 PetscInt M; 47 48 PetscFunctionBegin; 49 PetscCheck(!a->left && !a->right && irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented"); 50 if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat)); 51 PetscCall(MatGetSize(B, &M, NULL)); 52 PetscCall(PetscMalloc1(n, &row)); 53 PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0])); 54 PetscCall(ISSetIdentity(row[0])); 55 for (M = 1; M < n; ++M) row[M] = row[0]; 56 PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba)); 57 for (M = 0; M < n; ++M) { 58 PetscCall(MatCreateNormalHermitian(suba[M], *submat + M)); 59 ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale; 60 } 61 PetscCall(ISDestroy(&row[0])); 62 PetscCall(PetscFree(row)); 63 PetscCall(MatDestroySubMatrices(n, &suba)); 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B) { 68 Mat_Normal *a = (Mat_Normal *)A->data; 69 Mat C, Aa = a->A; 70 IS row; 71 72 PetscFunctionBegin; 73 PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same"); 74 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row)); 75 PetscCall(ISSetIdentity(row)); 76 PetscCall(MatPermute(Aa, row, colp, &C)); 77 PetscCall(ISDestroy(&row)); 78 PetscCall(MatCreateNormalHermitian(C, B)); 79 PetscCall(MatDestroy(&C)); 80 PetscFunctionReturn(0); 81 } 82 83 PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B) { 84 Mat_Normal *a = (Mat_Normal *)A->data; 85 Mat C; 86 87 PetscFunctionBegin; 88 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 89 PetscCall(MatDuplicate(a->A, op, &C)); 90 PetscCall(MatCreateNormalHermitian(C, B)); 91 PetscCall(MatDestroy(&C)); 92 if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale; 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str) { 97 Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data; 98 99 PetscFunctionBegin; 100 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 101 PetscCall(MatCopy(a->A, b->A, str)); 102 b->scale = a->scale; 103 PetscCall(VecDestroy(&b->left)); 104 PetscCall(VecDestroy(&b->right)); 105 PetscCall(VecDestroy(&b->leftwork)); 106 PetscCall(VecDestroy(&b->rightwork)); 107 PetscFunctionReturn(0); 108 } 109 110 PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y) { 111 Mat_Normal *Na = (Mat_Normal *)N->data; 112 Vec in; 113 114 PetscFunctionBegin; 115 in = x; 116 if (Na->right) { 117 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 118 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 119 in = Na->rightwork; 120 } 121 PetscCall(MatMult(Na->A, in, Na->w)); 122 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y)); 123 if (Na->left) PetscCall(VecPointwiseMult(y, Na->left, y)); 124 PetscCall(VecScale(y, Na->scale)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 129 Mat_Normal *Na = (Mat_Normal *)N->data; 130 Vec in; 131 132 PetscFunctionBegin; 133 in = v1; 134 if (Na->right) { 135 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 136 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 137 in = Na->rightwork; 138 } 139 PetscCall(MatMult(Na->A, in, Na->w)); 140 PetscCall(VecScale(Na->w, Na->scale)); 141 if (Na->left) { 142 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3)); 143 PetscCall(VecPointwiseMult(v3, Na->left, v3)); 144 PetscCall(VecAXPY(v3, 1.0, v2)); 145 } else { 146 PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3)); 147 } 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y) { 152 Mat_Normal *Na = (Mat_Normal *)N->data; 153 Vec in; 154 155 PetscFunctionBegin; 156 in = x; 157 if (Na->left) { 158 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 159 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 160 in = Na->leftwork; 161 } 162 PetscCall(MatMult(Na->A, in, Na->w)); 163 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y)); 164 if (Na->right) PetscCall(VecPointwiseMult(y, Na->right, y)); 165 PetscCall(VecScale(y, Na->scale)); 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 170 Mat_Normal *Na = (Mat_Normal *)N->data; 171 Vec in; 172 173 PetscFunctionBegin; 174 in = v1; 175 if (Na->left) { 176 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 177 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 178 in = Na->leftwork; 179 } 180 PetscCall(MatMult(Na->A, in, Na->w)); 181 PetscCall(VecScale(Na->w, Na->scale)); 182 if (Na->right) { 183 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3)); 184 PetscCall(VecPointwiseMult(v3, Na->right, v3)); 185 PetscCall(VecAXPY(v3, 1.0, v2)); 186 } else { 187 PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3)); 188 } 189 PetscFunctionReturn(0); 190 } 191 192 PetscErrorCode MatDestroy_NormalHermitian(Mat N) { 193 Mat_Normal *Na = (Mat_Normal *)N->data; 194 195 PetscFunctionBegin; 196 PetscCall(MatDestroy(&Na->A)); 197 PetscCall(MatDestroy(&Na->D)); 198 PetscCall(VecDestroy(&Na->w)); 199 PetscCall(VecDestroy(&Na->left)); 200 PetscCall(VecDestroy(&Na->right)); 201 PetscCall(VecDestroy(&Na->leftwork)); 202 PetscCall(VecDestroy(&Na->rightwork)); 203 PetscCall(PetscFree(N->data)); 204 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMatHermitian_C", NULL)); 205 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL)); 206 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL)); 207 PetscFunctionReturn(0); 208 } 209 210 /* 211 Slow, nonscalable version 212 */ 213 PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v) { 214 Mat_Normal *Na = (Mat_Normal *)N->data; 215 Mat A = Na->A; 216 PetscInt i, j, rstart, rend, nnz; 217 const PetscInt *cols; 218 PetscScalar *diag, *work, *values; 219 const PetscScalar *mvalues; 220 221 PetscFunctionBegin; 222 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work)); 223 PetscCall(PetscArrayzero(work, A->cmap->N)); 224 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 225 for (i = rstart; i < rend; i++) { 226 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues)); 227 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]); 228 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues)); 229 } 230 PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N))); 231 rstart = N->cmap->rstart; 232 rend = N->cmap->rend; 233 PetscCall(VecGetArray(v, &values)); 234 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart)); 235 PetscCall(VecRestoreArray(v, &values)); 236 PetscCall(PetscFree2(diag, work)); 237 PetscCall(VecScale(v, Na->scale)); 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D) { 242 Mat_Normal *Na = (Mat_Normal *)N->data; 243 Mat M, A = Na->A; 244 245 PetscFunctionBegin; 246 PetscCall(MatGetDiagonalBlock(A, &M)); 247 PetscCall(MatCreateNormalHermitian(M, &Na->D)); 248 *D = Na->D; 249 PetscFunctionReturn(0); 250 } 251 252 PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M) { 253 Mat_Normal *Aa = (Mat_Normal *)A->data; 254 255 PetscFunctionBegin; 256 *M = Aa->A; 257 PetscFunctionReturn(0); 258 } 259 260 /*@ 261 MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN` 262 263 Logically collective on A 264 265 Input Parameter: 266 . A - the `MATNORMALHERMITIAN` matrix 267 268 Output Parameter: 269 . M - the matrix object stored inside A 270 271 Level: intermediate 272 273 .seealso: `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()` 274 @*/ 275 PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M) { 276 PetscFunctionBegin; 277 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 278 PetscValidType(A, 1); 279 PetscValidPointer(M, 2); 280 PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M)); 281 PetscFunctionReturn(0); 282 } 283 284 PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 285 Mat_Normal *Aa = (Mat_Normal *)A->data; 286 Mat B, conjugate; 287 PetscInt m, n, M, N; 288 289 PetscFunctionBegin; 290 PetscCall(MatGetSize(A, &M, &N)); 291 PetscCall(MatGetLocalSize(A, &m, &n)); 292 if (reuse == MAT_REUSE_MATRIX) { 293 B = *newmat; 294 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B)); 295 } else { 296 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B)); 297 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 298 PetscCall(MatProductSetFromOptions(B)); 299 PetscCall(MatProductSymbolic(B)); 300 PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE)); 301 } 302 if (PetscDefined(USE_COMPLEX)) { 303 PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate)); 304 PetscCall(MatConjugate(conjugate)); 305 PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B)); 306 } 307 PetscCall(MatProductNumeric(B)); 308 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate)); 309 if (reuse == MAT_INPLACE_MATRIX) { 310 PetscCall(MatHeaderReplace(A, &B)); 311 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 312 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat)); 313 PetscFunctionReturn(0); 314 } 315 316 /*@ 317 MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A. 318 319 Collective on A 320 321 Input Parameter: 322 . A - the (possibly rectangular complex) matrix 323 324 Output Parameter: 325 . N - the matrix that represents (A*)'*A 326 327 Level: intermediate 328 329 Note: 330 The product (A*)'*A is NOT actually formed! Rather the new matrix 331 object performs the matrix-vector product, `MatMult()`, by first multiplying by 332 A and then (A*)' 333 334 .seealso: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()` 335 @*/ 336 PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N) { 337 PetscInt m, n; 338 Mat_Normal *Na; 339 VecType vtype; 340 341 PetscFunctionBegin; 342 PetscCall(MatGetLocalSize(A, &m, &n)); 343 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 344 PetscCall(MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE)); 345 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN)); 346 PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap)); 347 PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap)); 348 349 PetscCall(PetscNew(&Na)); 350 (*N)->data = (void *)Na; 351 PetscCall(PetscObjectReference((PetscObject)A)); 352 Na->A = A; 353 Na->scale = 1.0; 354 355 PetscCall(MatCreateVecs(A, NULL, &Na->w)); 356 357 (*N)->ops->destroy = MatDestroy_NormalHermitian; 358 (*N)->ops->mult = MatMult_NormalHermitian; 359 (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal; 360 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal; 361 (*N)->ops->multadd = MatMultHermitianAdd_Normal; 362 (*N)->ops->getdiagonal = MatGetDiagonal_NormalHermitian; 363 (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_NormalHermitian; 364 (*N)->ops->scale = MatScale_NormalHermitian; 365 (*N)->ops->diagonalscale = MatDiagonalScale_NormalHermitian; 366 (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian; 367 (*N)->ops->permute = MatPermute_NormalHermitian; 368 (*N)->ops->duplicate = MatDuplicate_NormalHermitian; 369 (*N)->ops->copy = MatCopy_NormalHermitian; 370 (*N)->assembled = PETSC_TRUE; 371 (*N)->preallocated = PETSC_TRUE; 372 373 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMatHermitian_C", MatNormalGetMat_NormalHermitian)); 374 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ)); 375 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ)); 376 PetscCall(MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE)); 377 PetscCall(MatGetVecType(A, &vtype)); 378 PetscCall(MatSetVecType(*N, vtype)); 379 #if defined(PETSC_HAVE_DEVICE) 380 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 381 #endif 382 PetscFunctionReturn(0); 383 } 384