Actual source code: stshellmat.c
slepc-3.21.1 2024-04-26
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: */
10: /*
11: Subroutines that implement various operations of the matrix associated with
12: the shift-and-invert technique for eigenvalue problems
13: */
15: #include <slepc/private/stimpl.h>
17: typedef struct {
18: PetscScalar alpha;
19: PetscScalar *coeffs;
20: ST st;
21: Vec z;
22: PetscInt nmat;
23: PetscInt *matIdx;
24: } ST_MATSHELL;
26: PetscErrorCode STMatShellShift(Mat A,PetscScalar alpha)
27: {
28: ST_MATSHELL *ctx;
30: PetscFunctionBegin;
31: PetscCall(MatShellGetContext(A,&ctx));
32: ctx->alpha = alpha;
33: PetscCall(PetscObjectStateIncrease((PetscObject)A));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: /*
38: For i=0:nmat-1 computes y = (sum_i (coeffs[i]*alpha^i*st->A[idx[i]]))x
39: If null coeffs computes with coeffs[i]=1.0
40: */
41: static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
42: {
43: ST_MATSHELL *ctx;
44: ST st;
45: PetscInt i;
46: PetscScalar t=1.0,c;
48: PetscFunctionBegin;
49: PetscCall(MatShellGetContext(A,&ctx));
50: st = ctx->st;
51: PetscCall(MatMult(st->A[ctx->matIdx[0]],x,y));
52: if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
53: if (ctx->alpha!=0.0) {
54: for (i=1;i<ctx->nmat;i++) {
55: PetscCall(MatMult(st->A[ctx->matIdx[i]],x,ctx->z));
56: t *= ctx->alpha;
57: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
58: PetscCall(VecAXPY(y,c,ctx->z));
59: }
60: if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
66: {
67: ST_MATSHELL *ctx;
68: ST st;
69: PetscInt i;
70: PetscScalar t=1.0,c;
72: PetscFunctionBegin;
73: PetscCall(MatShellGetContext(A,&ctx));
74: st = ctx->st;
75: PetscCall(MatMultTranspose(st->A[ctx->matIdx[0]],x,y));
76: if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
77: if (ctx->alpha!=0.0) {
78: for (i=1;i<ctx->nmat;i++) {
79: PetscCall(MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z));
80: t *= ctx->alpha;
81: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
82: PetscCall(VecAXPY(y,c,ctx->z));
83: }
84: if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
85: }
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: #if defined(PETSC_USE_COMPLEX)
90: static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A,Vec x,Vec y)
91: {
92: ST_MATSHELL *ctx;
93: ST st;
94: PetscInt i;
95: PetscScalar t=1.0,c;
97: PetscFunctionBegin;
98: PetscCall(MatShellGetContext(A,&ctx));
99: st = ctx->st;
100: PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[0]],x,y));
101: if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,PetscConj(ctx->coeffs[0])));
102: if (ctx->alpha!=0.0) {
103: for (i=1;i<ctx->nmat;i++) {
104: PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[i]],x,ctx->z));
105: t *= ctx->alpha;
106: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
107: PetscCall(VecAXPY(y,PetscConj(c),ctx->z));
108: }
109: if (ctx->nmat==1) PetscCall(VecAXPY(y,PetscConj(ctx->alpha),x)); /* y = (A + alpha*I) x */
110: }
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
113: #endif
115: static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
116: {
117: ST_MATSHELL *ctx;
118: ST st;
119: Vec diagb;
120: PetscInt i;
121: PetscScalar t=1.0,c;
123: PetscFunctionBegin;
124: PetscCall(MatShellGetContext(A,&ctx));
125: st = ctx->st;
126: PetscCall(MatGetDiagonal(st->A[ctx->matIdx[0]],diag));
127: if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(diag,ctx->coeffs[0]));
128: if (ctx->alpha!=0.0) {
129: if (ctx->nmat==1) PetscCall(VecShift(diag,ctx->alpha)); /* y = (A + alpha*I) x */
130: else {
131: PetscCall(VecDuplicate(diag,&diagb));
132: for (i=1;i<ctx->nmat;i++) {
133: PetscCall(MatGetDiagonal(st->A[ctx->matIdx[i]],diagb));
134: t *= ctx->alpha;
135: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
136: PetscCall(VecAYPX(diag,c,diagb));
137: }
138: PetscCall(VecDestroy(&diagb));
139: }
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode MatDestroy_Shell(Mat A)
145: {
146: ST_MATSHELL *ctx;
148: PetscFunctionBegin;
149: PetscCall(MatShellGetContext(A,&ctx));
150: PetscCall(VecDestroy(&ctx->z));
151: PetscCall(PetscFree(ctx->matIdx));
152: PetscCall(PetscFree(ctx->coeffs));
153: PetscCall(PetscFree(ctx));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
158: {
159: PetscInt n,m,N,M,i;
160: PetscBool has=PETSC_FALSE,hasA,hasB;
161: ST_MATSHELL *ctx;
163: PetscFunctionBegin;
164: PetscCall(MatGetSize(st->A[0],&M,&N));
165: PetscCall(MatGetLocalSize(st->A[0],&m,&n));
166: PetscCall(PetscNew(&ctx));
167: ctx->st = st;
168: ctx->alpha = alpha;
169: ctx->nmat = matIdx?nmat:st->nmat;
170: PetscCall(PetscMalloc1(ctx->nmat,&ctx->matIdx));
171: if (matIdx) {
172: for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
173: } else {
174: ctx->matIdx[0] = 0;
175: if (ctx->nmat>1) ctx->matIdx[1] = 1;
176: }
177: if (coeffs) {
178: PetscCall(PetscMalloc1(ctx->nmat,&ctx->coeffs));
179: for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
180: }
181: PetscCall(MatCreateVecs(st->A[0],&ctx->z,NULL));
182: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat));
183: PetscCall(MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell));
184: PetscCall(MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
185: #if defined(PETSC_USE_COMPLEX)
186: PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_Shell));
187: #else
188: PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
189: #endif
190: PetscCall(MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell));
192: PetscCall(MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA));
193: if (st->nmat>1) {
194: has = hasA;
195: for (i=1;i<ctx->nmat;i++) {
196: PetscCall(MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB));
197: has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
198: }
199: }
200: if ((hasA && st->nmat==1) || has) PetscCall(MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }