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