Actual source code: matstruct.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/slepcimpl.h>
13: /*@
14: MatCreateBSE - Create a matrix that can be used to define a structured eigenvalue
15: problem of type BSE (Bethe-Salpeter Equation).
17: Collective
19: Input Parameters:
20: + R - matrix for the diagonal block (resonant)
21: - C - matrix for the off-diagonal block (coupling)
23: Output Parameter:
24: . H - the resulting matrix
26: Notes:
27: The resulting matrix has the block form $$H = \begin{bmatrix} R & C \\ -C^* & -R^T \end{bmatrix},$$
28: where $R$ is assumed
29: to be (complex) Hermitian and $C$ complex symmetric. Note that this function does
30: not check these properties, so if the matrices provided by the user do not satisfy
31: them, then the solver will not behave as expected.
33: The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
34: `EPSSetOperators()` for the case that the problem type is `EPS_BSE`. Note that the user
35: cannot just build a matrix with the required structure, it must be done via this
36: function.
38: In the current implementation, `H` is a `MATNEST` matrix, where `R` and `C` form the top
39: block row, while the bottom block row is composed of matrices of type
40: `MATTRANSPOSEVIRTUAL` and `MATHERMITIANTRANSPOSEVIRTUAL` scaled by -1.
42: Level: intermediate
44: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateHamiltonian()`, `MatCreateLREP()`
45: @*/
46: PetscErrorCode MatCreateBSE(Mat R,Mat C,Mat *H)
47: {
48: PetscInt Mr,Mc,Nr,Nc,mr,mc,nr,nc;
49: Mat block[4] = { R, C, NULL, NULL };
50: SlepcMatStruct mctx;
52: PetscFunctionBegin;
55: PetscCheckSameTypeAndComm(R,1,C,2);
56: PetscAssertPointer(H,3);
58: /* check sizes */
59: PetscCall(MatGetSize(R,&Mr,&Nr));
60: PetscCall(MatGetLocalSize(R,&mr,&nr));
61: PetscCall(MatGetSize(C,&Mc,&Nc));
62: PetscCall(MatGetLocalSize(C,&mc,&nc));
63: PetscCheck(Mc==Mr && mc==mr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
64: PetscCheck(Nc==Nr && nc==nr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");
66: /* bottom block row */
67: PetscCall(MatCreateHermitianTranspose(C,&block[2]));
68: PetscCall(MatScale(block[2],-1.0));
69: PetscCall(MatCreateTranspose(R,&block[3]));
70: PetscCall(MatScale(block[3],-1.0));
72: /* create nest matrix and compose context */
73: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R),2,NULL,2,NULL,block,H));
74: PetscCall(MatDestroy(&block[2]));
75: PetscCall(MatDestroy(&block[3]));
77: PetscCall(PetscNew(&mctx));
78: mctx->cookie = SLEPC_MAT_STRUCT_BSE;
79: PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mctx,PetscCtxDestroyDefault));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: /*@
84: MatCreateHamiltonian - Create a matrix that can be used to define a structured
85: eigenvalue problem of Hamiltonian type.
87: Collective
89: Input Parameters:
90: + A - matrix for the (0,0) block
91: . B - matrix for the (0,1) block, must be real symmetric or complex Hermitian
92: - C - matrix for the (1,0) block, must be real symmetric or complex Hermitian
94: Output Parameter:
95: . H - the resulting matrix
97: Notes:
98: The resulting matrix has the block form $$H = \begin{bmatrix} A & B \\ C & -A^* \end{bmatrix},$$
99: where $B$ and $C$ are
100: assumed to be symmetric in the real case or Hermitian in the complex case. Note
101: that this function does not check this property, so if the matrices provided by
102: the user do not satisfy it, then the solver will not behave as expected.
104: The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
105: `EPSSetOperators()` for the case that the problem type is `EPS_HAMILT`. Note that the
106: user cannot just build a matrix with the required structure, it must be done via
107: this function.
109: In the current implementation, `H` is a `MATNEST` matrix, where the (1,1) block is
110: a matrix of type `MATHERMITIANTRANSPOSEVIRTUAL` obtained from `A` and scaled by -1.
112: Level: intermediate
114: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateBSE()`, `MatCreateLREP()`
115: @*/
116: PetscErrorCode MatCreateHamiltonian(Mat A,Mat B,Mat C,Mat *H)
117: {
118: PetscInt Ma,Mb,Mc,Na,Nb,Nc,ma,mb,mc,na,nb,nc;
119: Mat block[4] = { A, B, C, NULL };
120: SlepcMatStruct mctx;
122: PetscFunctionBegin;
126: PetscCheckSameTypeAndComm(A,1,B,2);
127: PetscCheckSameTypeAndComm(A,1,C,3);
128: PetscAssertPointer(H,4);
130: /* check sizes */
131: PetscCall(MatGetSize(A,&Ma,&Na));
132: PetscCall(MatGetLocalSize(A,&ma,&na));
133: PetscCall(MatGetSize(B,&Mb,&Nb));
134: PetscCall(MatGetLocalSize(B,&mb,&nb));
135: PetscCall(MatGetSize(C,&Mc,&Nc));
136: PetscCall(MatGetLocalSize(C,&mc,&nc));
137: PetscCheck(Mb==Ma && mb==ma,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
138: PetscCheck(Nc==Na && nc==na,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");
140: /* (1,1) block */
141: PetscCall(MatCreateHermitianTranspose(A,&block[3]));
142: PetscCall(MatScale(block[3],-1.0));
144: /* create nest matrix and compose context */
145: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A),2,NULL,2,NULL,block,H));
146: PetscCall(MatDestroy(&block[3]));
148: PetscCall(PetscNew(&mctx));
149: mctx->cookie = SLEPC_MAT_STRUCT_HAMILT;
150: PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mctx,PetscCtxDestroyDefault));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*@
155: MatCreateLREP - Create a matrix that can be used to define a structured Linear
156: Response eigenvalue problem.
158: Collective
160: Input Parameters:
161: + AK - matrix for the diagonal block ($A$) or the top block ($K$)
162: . BM - matrix for the off-diagonal block ($B$) or the bottom block ($M$)
163: - red - whether the reduced form should be built
165: Output Parameter:
166: . H - the resulting matrix
168: Notes:
169: The resulting matrix has the block form $$H = \begin{bmatrix} A & B \\ -B & -A \end{bmatrix},$$
170: where both $A$ and $B$ are assumed to be real symmetric. An alternative form,
171: called reduced form (`red`=`PETSC_TRUE`) is $$H = \begin{bmatrix} 0 & K \\ M & 0 \end{bmatrix},$$
172: where both $K$ and $M$ are also assumed to be real symmetric. Note that this function
173: does not check these properties, so if the matrices provided by the user do not
174: satisfy them, then the solver will not behave as expected.
176: The reduced form can be obtained from the original formulation by setting $K=A-B$
177: and $M=A+B$. However, $K$ and $M$ need not have such relation, in which case the
178: computed eigenvalues will be $\pm\lambda_i$, where $\lambda_i^2$ are the eigenvalues
179: of the _product eigenvalue problem_ $KM$ (or $MK$).
181: The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
182: `EPSSetOperators()` for the case that the problem type is `EPS_LREP`. Note that the
183: user cannot just build a matrix with the required structure, it must be done via
184: this function.
186: In the current implementation, `H` is a `MATNEST` matrix, where `A` and `B` form the
187: top block row, while the bottom block row is composed of shell matrices scaled by -1.
188: Similarly, in the reduced formulation, the `MATNEST` has two zero blocks.
190: Internally, the eigensolver will work with the reduced form, and it will carry out
191: the conversion transparently if the first form is provided by the user.
193: Level: intermediate
195: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateBSE()`, `MatCreateHamiltonian()`
196: @*/
197: PetscErrorCode MatCreateLREP(Mat AK,Mat BM,PetscBool red,Mat *H)
198: {
199: PetscInt Ma,Mb,Na,Nb,ma,mb,na,nb;
200: SlepcMatStruct mbtx;
202: PetscFunctionBegin;
205: PetscCheckSameTypeAndComm(AK,1,BM,2);
207: PetscAssertPointer(H,4);
209: /* check sizes */
210: PetscCall(MatGetSize(AK,&Ma,&Na));
211: PetscCall(MatGetLocalSize(AK,&ma,&na));
212: PetscCall(MatGetSize(BM,&Mb,&Nb));
213: PetscCall(MatGetLocalSize(BM,&mb,&nb));
214: PetscCheck(Mb==Ma && mb==ma,PetscObjectComm((PetscObject)AK),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
215: PetscCheck(Nb==Na && nb==na,PetscObjectComm((PetscObject)AK),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");
217: if (red) { /* reduced form */
219: Mat block[4] = { NULL, AK, BM, NULL };
221: /* create nest matrix */
222: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)AK),2,NULL,2,NULL,block,H));
224: } else { /* original LREP form */
226: Mat block[4] = { AK, BM, NULL, NULL };
227: const PetscScalar scal[] = { -1.0 };
229: /* bottom block row, using MATCOMPOSITE */
230: PetscCall(MatCreate(PetscObjectComm((PetscObject)AK),&block[2]));
231: PetscCall(MatSetSizes(block[2],mb,nb,Mb,Nb));
232: PetscCall(MatSetType(block[2],MATCOMPOSITE));
233: PetscCall(MatCompositeAddMat(block[2],BM));
234: PetscCall(MatAssemblyBegin(block[2],MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(block[2],MAT_FINAL_ASSEMBLY));
236: PetscCall(MatCompositeSetScalings(block[2],scal));
238: PetscCall(MatCreate(PetscObjectComm((PetscObject)AK),&block[3]));
239: PetscCall(MatSetSizes(block[3],ma,na,Ma,Na));
240: PetscCall(MatSetType(block[3],MATCOMPOSITE));
241: PetscCall(MatCompositeAddMat(block[3],AK));
242: PetscCall(MatAssemblyBegin(block[3],MAT_FINAL_ASSEMBLY));
243: PetscCall(MatAssemblyEnd(block[3],MAT_FINAL_ASSEMBLY));
244: PetscCall(MatCompositeSetScalings(block[3],scal));
246: /* create nest matrix */
247: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)AK),2,NULL,2,NULL,block,H));
248: PetscCall(MatDestroy(&block[2]));
249: PetscCall(MatDestroy(&block[3]));
251: }
253: /* compose context */
254: PetscCall(PetscNew(&mbtx));
255: mbtx->cookie = SLEPC_MAT_STRUCT_LREP;
256: PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mbtx,PetscCtxDestroyDefault));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }