1: #include <../src/mat/impls/nest/matnestimpl.h>

  3: static PetscErrorCode MatMultiShiftDestroy(PetscCtxRt ctx)
  4: {
  5:   Mat_MultiShift mctx = *(Mat_MultiShift *)ctx;

  7:   PetscFunctionBegin;
  8:   PetscCall(MatDestroy(&mctx->K));
  9:   PetscCall(MatDestroy(&mctx->M));
 10:   PetscCall(PetscFree(mctx->sigma));
 11:   PetscCall(PetscFree(mctx->cmplx));
 12:   PetscCall(PetscFree(mctx));
 13:   PetscFunctionReturn(PETSC_SUCCESS);
 14: }

 16: /* Build G = K + sigma M, explicitly or not */
 17: PetscErrorCode MatMultiShiftBuildShiftedMatrix_Internal(Mat K, PetscScalar sigma, Mat M, MatStructure str, PetscBool explicitmat, Mat *G)
 18: {
 19:   PetscScalar scal[] = {1.0, 1.0};

 21:   PetscFunctionBegin;
 22:   if (explicitmat) {
 23:     PetscCall(MatDuplicate(K, MAT_COPY_VALUES, G));
 24:     if (M) PetscCall(MatAXPY(*G, sigma, M, str));
 25:     else PetscCall(MatShift(*G, sigma));
 26:   } else {
 27:     PetscInt Mk, Nk, mk, nk;

 29:     PetscCall(MatGetSize(K, &Mk, &Nk));
 30:     PetscCall(MatGetLocalSize(K, &mk, &nk));
 31:     PetscCall(MatCreate(PetscObjectComm((PetscObject)K), G));
 32:     PetscCall(MatSetSizes(*G, mk, nk, Mk, Nk));
 33:     PetscCall(MatSetType(*G, MATCOMPOSITE));
 34:     PetscCall(MatCompositeAddMat(*G, K));
 35:     if (M) PetscCall(MatCompositeAddMat(*G, M));
 36:     PetscCall(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY));
 37:     PetscCall(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY));
 38:     if (M) {
 39:       scal[1] = sigma;
 40:       PetscCall(MatCompositeSetScalings(*G, scal));
 41:     } else PetscCall(MatShift(*G, sigma));
 42:   }
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: #if !PetscDefined(USE_COMPLEX)
 47: /* Build G = sigma M, explicitly or not, where M may be I; K is used for dimensions only */
 48: static PetscErrorCode BuildScaledMatrix(Mat K, PetscScalar sigma, Mat M, PetscBool explicitmat, Mat *G)
 49: {
 50:   PetscFunctionBegin;
 51:   if (explicitmat) {
 52:     if (M) {
 53:       PetscCall(MatDuplicate(M, MAT_COPY_VALUES, G));
 54:       PetscCall(MatScale(*G, sigma));
 55:     } else { // M=I
 56:       PetscCall(MatDuplicate(K, MAT_DO_NOT_COPY_VALUES, G));
 57:       PetscCall(MatZeroEntries(*G));
 58:       PetscCall(MatShift(*G, sigma));
 59:     }
 60:   } else {
 61:     PetscInt Mk, Nk, mk, nk;

 63:     PetscCall(MatGetSize(K, &Mk, &Nk));
 64:     PetscCall(MatGetLocalSize(K, &mk, &nk));
 65:     if (M) {
 66:       PetscCall(MatCreate(PetscObjectComm((PetscObject)K), G));
 67:       PetscCall(MatSetSizes(*G, mk, nk, Mk, Nk));
 68:       PetscCall(MatSetType(*G, MATCOMPOSITE));
 69:       PetscCall(MatCompositeAddMat(*G, M));
 70:       PetscCall(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY));
 71:       PetscCall(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY));
 72:       PetscCall(MatCompositeSetScalings(*G, &sigma));
 73:     } else PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)K), mk, nk, Mk, Nk, sigma, G));
 74:   }
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }
 77: #endif

 79: /*@
 80:   MatCreateNestFromMultipleShifts - Creates a `MATNEST` matrix that represents a family of shifted matrices
 81:   $K + \sigma_i M$ for a number of shifts $\sigma_i$.

 83:   Collective

 85:   Input Parameters:
 86: + K               - the first `Mat` (stiffness) forming the shifted matrices
 87: . nshift          - number of shifts
 88: . sigma           - array of shifts $\sigma_i$, its length is `nshift`
 89: . sigma_imaginary - imaginary parts of the shifts $\sigma_i$ in case of complex-conjugate pairs (can be `NULL`); only used when `PetscScalar` is `PetscReal`
 90: . M               - the second `Mat` (mass) forming the shifted matrices (if `NULL` the identity matrix is assumed)
 91: . explicitmat     - whether the shifted matrices should be built explicitly or not
 92: - str             - `MatStructure` flag

 94:   Output Parameter:
 95: . A - the resulting matrix

 97:   Level: intermediate

 99:   Notes:
100:   This is intended for solving a family of shifted linear systems, $(K + \sigma_i M) x_i = b$,
101:   where in some applications $M = I$. This function returns a `MATNEST` `A` whose `nshift`
102:   diagonal blocks are the matrices $K + \sigma_i M$, either built explicitly or not, depending
103:   on the `explicitmat` argument. If not explicit, the diagonal blocks are created as `MATCOMPOSITE`.
104:   If `explicitmat` is true, and $M$ is not the identity, then $K + \sigma_i M$ is built with a
105:   call to `MatAXPY()`, where the flag `str` is used to indicate the relation between the sparsity
106:   patterns of $K$ and $M$. The flag is also stored in `A` and reused for the same purpose by
107:   solvers that build a shifted matrix explicitly from `A`, such as `KSPEKSM`, so it must describe
108:   the relation between the two patterns even when `explicitmat` is `PETSC_FALSE`.

110:   To solve all the shifted linear systems simultaneously, pass this matrix to a
111:   `KSP` solver such as `KSPEKSM`, along with compatible solution and right-hand side vectors. Since
112:   the right-hand side $b$ is the same for all the shifted linear systems, one can call
113:   `MatCreateVecNestFromMultipleShifts()` to easily create a nested `Vec` containing `nshift` references to $b$.

115:   When `PetscScalar` is `PetscReal` it is possible to provide complex conjugate pairs of shifts by
116:   passing a nonzero imaginary part in `sigma_imaginary[i]`. In that case `A` is no longer block
117:   diagonal but contains 2x2 diagonal blocks for each complex-conjugate pair.

119: .seealso: [](ch_matrices), `MatCreateVecNestFromMultipleShifts()`, `MatAXPY()`, `MATCOMPOSITE`, `KSP`, `KSPEKSM`
120: @*/
121: PetscErrorCode MatCreateNestFromMultipleShifts(Mat K, PetscInt nshift, const PetscScalar sigma[], const PetscScalar sigma_imaginary[], Mat M, PetscBool explicitmat, MatStructure str, Mat *A)
122: {
123:   PetscInt       i;
124:   Mat           *mats;
125:   Mat_MultiShift mctx;

127:   PetscFunctionBegin;
130:   PetscAssertPointer(sigma, 3);
131:   if (sigma_imaginary) PetscAssertPointer(sigma_imaginary, 4);
135:   PetscAssertPointer(A, 8);

137:   PetscCheck(nshift > 0, PetscObjectComm((PetscObject)K), PETSC_ERR_ARG_OUTOFRANGE, "The value nshift must be > 0");

139:   /* build context and array of shifted matrices */
140:   PetscCall(PetscNew(&mctx));
141:   PetscCall(PetscObjectReference((PetscObject)K));
142:   mctx->K = K;
143:   PetscCall(PetscObjectReference((PetscObject)M));
144:   mctx->M      = M;
145:   mctx->nshift = nshift;
146:   mctx->str    = str;
147:   PetscCall(PetscMalloc1(nshift, &mctx->sigma));
148:   if (!PetscDefined(USE_COMPLEX)) PetscCall(PetscCalloc1(nshift, &mctx->cmplx));
149:   PetscCall(PetscCalloc1(nshift * nshift, &mats));

151:   for (i = 0; i < nshift; i++) {
152: #if PetscDefined(USE_COMPLEX)
153:     mctx->sigma[i] = sigma[i];
154:     PetscCall(MatMultiShiftBuildShiftedMatrix_Internal(K, sigma[i], M, str, explicitmat, mats + i + i * nshift));
155: #else
156:     mctx->sigma[i] = sigma[i];
157:     PetscCall(MatMultiShiftBuildShiftedMatrix_Internal(K, mctx->sigma[i], M, str, explicitmat, mats + i + i * nshift));
158:     if (sigma_imaginary && sigma_imaginary[i] != 0.0) {
159:       PetscCheck(i < nshift - 1, PetscObjectComm((PetscObject)K), PETSC_ERR_ARG_WRONG, "The last shift of the array is complex; shifts must come in complex-conjugate pairs");
160:       PetscCheck(sigma[i + 1] == sigma[i] && sigma_imaginary[i + 1] == -sigma_imaginary[i], PetscObjectComm((PetscObject)K), PETSC_ERR_ARG_WRONG, "The shifts must either be real or form a (consecutive) complex-conjugate pair");
161:       mctx->sigma[i + 1] = sigma_imaginary[i];
162:       mctx->cmplx[i]     = PETSC_TRUE;
163:       // 2x2 block in the nested matrix
164:       PetscCall(PetscObjectReference((PetscObject)mats[i + i * nshift]));
165:       mats[i + 1 + (i + 1) * nshift] = mats[i + i * nshift];
166:       PetscCall(BuildScaledMatrix(K, sigma_imaginary[i], M, explicitmat, mats + i + (i + 1) * nshift));
167:       PetscCall(BuildScaledMatrix(K, -sigma_imaginary[i], M, explicitmat, mats + i + 1 + i * nshift));
168:       i++; // skip next shift
169:     }
170: #endif
171:   }

173:   /* build MATNEST */
174:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)K), nshift, NULL, nshift, NULL, mats, A));
175:   for (i = 0; i < nshift * nshift; i++) PetscCall(MatDestroy(&mats[i]));
176:   PetscCall(PetscFree(mats));
177:   PetscCall(MatNestSetVecType(*A, VECNEST));

179:   /* compose context */
180:   PetscCall(MatGetState(*A, &mctx->state)); // used to detect later changes such as MatScale()
181:   PetscCall(PetscObjectContainerCompose((PetscObject)*A, "MatMultiShift", mctx, MatMultiShiftDestroy));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: /*@
186:   MatCreateVecNestFromMultipleShifts - Creates a `VECNEST` that is compatible with a matrix
187:   created with `MatCreateNestFromMultipleShifts()`.

189:   Collective

191:   Input Parameters:
192: + A - a `Mat` created with `MatCreateNestFromMultipleShifts()`
193: - v - an optional vector (set to `NULL` if not needed)

195:   Output Parameter:
196: . vout - the resulting vector

198:   Level: intermediate

200:   Notes:
201:   The result is a `VECNEST` compatible with `A`, so that it can, e.g., be multiplied against.
202:   If the input vector `v` is passed, `vout` will contain `nshift` references to `v`, where
203:   `nshift` is the number of subvectors of `vout`. Hence, it is intended as a read-only right-hand
204:   side, not a solution vector.

206:   In real scalars, in the case of a complex-conjugate pair only the first block of the pair is
207:   set while the second is left zero. The reason is that in this case the second block represents
208:   the imaginary part, which is zero since `v` is real.

210: .seealso: [](ch_matrices), `MatCreateNestFromMultipleShifts()`
211: @*/
212: PetscErrorCode MatCreateVecNestFromMultipleShifts(Mat A, Vec v, Vec *vout)
213: {
214:   PetscInt       i;
215:   Mat_MultiShift mctx;

217:   PetscFunctionBegin;
220:   PetscAssertPointer(vout, 3);
221:   MatCheckMultiShift(A, &mctx);

223:   PetscCall(MatCreateVecs(A, NULL, vout));
224:   if (v)
225:     for (i = 0; i < mctx->nshift; i++) {
226:       PetscCall(VecNestSetSubVec(*vout, i, v));
227:       if (!PetscDefined(USE_COMPLEX) && mctx->cmplx[i]) i++; // complex shift, leave a zero block since b is real
228:     }
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }