Actual source code: stshellmat.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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_SHELLMAT;
26: PetscErrorCode STMatShellShift(Mat A,PetscScalar alpha)
27: {
29: ST_SHELLMAT *ctx;
32: MatShellGetContext(A,(void**)&ctx);
33: ctx->alpha = alpha;
34: return(0);
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: {
44: ST_SHELLMAT *ctx;
45: ST st;
46: PetscInt i;
47: PetscScalar t=1.0,c;
50: MatShellGetContext(A,(void**)&ctx);
51: st = ctx->st;
52: MatMult(st->A[ctx->matIdx[0]],x,y);
53: if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
54: VecScale(y,ctx->coeffs[0]);
55: }
56: if (ctx->alpha!=0.0) {
57: for (i=1;i<ctx->nmat;i++) {
58: MatMult(st->A[ctx->matIdx[i]],x,ctx->z);
59: t *= ctx->alpha;
60: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
61: VecAXPY(y,c,ctx->z);
62: }
63: if (ctx->nmat==1) { /* y = (A + alpha*I) x */
64: VecAXPY(y,ctx->alpha,x);
65: }
66: }
67: return(0);
68: }
70: static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
71: {
73: ST_SHELLMAT *ctx;
74: ST st;
75: PetscInt i;
76: PetscScalar t=1.0,c;
79: MatShellGetContext(A,(void**)&ctx);
80: st = ctx->st;
81: MatMultTranspose(st->A[ctx->matIdx[0]],x,y);
82: if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
83: VecScale(y,ctx->coeffs[0]);
84: }
85: if (ctx->alpha!=0.0) {
86: for (i=1;i<ctx->nmat;i++) {
87: MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z);
88: t *= ctx->alpha;
89: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
90: VecAXPY(y,c,ctx->z);
91: }
92: if (ctx->nmat==1) { /* y = (A + alpha*I) x */
93: VecAXPY(y,ctx->alpha,x);
94: }
95: }
96: return(0);
97: }
99: static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
100: {
102: ST_SHELLMAT *ctx;
103: ST st;
104: Vec diagb;
105: PetscInt i;
106: PetscScalar t=1.0,c;
109: MatShellGetContext(A,(void**)&ctx);
110: st = ctx->st;
111: MatGetDiagonal(st->A[ctx->matIdx[0]],diag);
112: if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
113: VecScale(diag,ctx->coeffs[0]);
114: }
115: if (ctx->alpha!=0.0) {
116: if (ctx->nmat==1) { /* y = (A + alpha*I) x */
117: VecShift(diag,ctx->alpha);
118: } else {
119: VecDuplicate(diag,&diagb);
120: for (i=1;i<ctx->nmat;i++) {
121: MatGetDiagonal(st->A[ctx->matIdx[i]],diagb);
122: t *= ctx->alpha;
123: c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
124: VecAYPX(diag,c,diagb);
125: }
126: VecDestroy(&diagb);
127: }
128: }
129: return(0);
130: }
132: static PetscErrorCode MatDestroy_Shell(Mat A)
133: {
135: ST_SHELLMAT *ctx;
138: MatShellGetContext(A,(void**)&ctx);
139: VecDestroy(&ctx->z);
140: PetscFree(ctx->matIdx);
141: PetscFree(ctx->coeffs);
142: PetscFree(ctx);
143: return(0);
144: }
146: PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
147: {
149: PetscInt n,m,N,M,i;
150: PetscBool has=PETSC_FALSE,hasA,hasB;
151: ST_SHELLMAT *ctx;
154: MatGetSize(st->A[0],&M,&N);
155: MatGetLocalSize(st->A[0],&m,&n);
156: PetscNew(&ctx);
157: ctx->st = st;
158: ctx->alpha = alpha;
159: ctx->nmat = matIdx?nmat:st->nmat;
160: PetscMalloc1(ctx->nmat,&ctx->matIdx);
161: if (matIdx) {
162: for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
163: } else {
164: ctx->matIdx[0] = 0;
165: if (ctx->nmat>1) ctx->matIdx[1] = 1;
166: }
167: if (coeffs) {
168: PetscMalloc1(ctx->nmat,&ctx->coeffs);
169: for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
170: }
171: MatCreateVecs(st->A[0],&ctx->z,NULL);
172: MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat);
173: MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell);
174: MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell);
175: MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell);
177: MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA);
178: if (st->nmat>1) {
179: has = hasA;
180: for (i=1;i<ctx->nmat;i++) {
181: MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB);
182: has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
183: }
184: }
185: if ((hasA && st->nmat==1) || has) {
186: MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);
187: }
188: return(0);
189: }