Actual source code: stshellmat.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  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: static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
 90: {
 91:   ST_MATSHELL    *ctx;
 92:   ST             st;
 93:   Vec            diagb;
 94:   PetscInt       i;
 95:   PetscScalar    t=1.0,c;

 97:   PetscFunctionBegin;
 98:   PetscCall(MatShellGetContext(A,&ctx));
 99:   st = ctx->st;
100:   PetscCall(MatGetDiagonal(st->A[ctx->matIdx[0]],diag));
101:   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(diag,ctx->coeffs[0]));
102:   if (ctx->alpha!=0.0) {
103:     if (ctx->nmat==1) PetscCall(VecShift(diag,ctx->alpha)); /* y = (A + alpha*I) x */
104:     else {
105:       PetscCall(VecDuplicate(diag,&diagb));
106:       for (i=1;i<ctx->nmat;i++) {
107:         PetscCall(MatGetDiagonal(st->A[ctx->matIdx[i]],diagb));
108:         t *= ctx->alpha;
109:         c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
110:         PetscCall(VecAYPX(diag,c,diagb));
111:       }
112:       PetscCall(VecDestroy(&diagb));
113:     }
114:   }
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode MatDestroy_Shell(Mat A)
119: {
120:   ST_MATSHELL    *ctx;

122:   PetscFunctionBegin;
123:   PetscCall(MatShellGetContext(A,&ctx));
124:   PetscCall(VecDestroy(&ctx->z));
125:   PetscCall(PetscFree(ctx->matIdx));
126:   PetscCall(PetscFree(ctx->coeffs));
127:   PetscCall(PetscFree(ctx));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
132: {
133:   PetscInt       n,m,N,M,i;
134:   PetscBool      has=PETSC_FALSE,hasA,hasB;
135:   ST_MATSHELL    *ctx;

137:   PetscFunctionBegin;
138:   PetscCall(MatGetSize(st->A[0],&M,&N));
139:   PetscCall(MatGetLocalSize(st->A[0],&m,&n));
140:   PetscCall(PetscNew(&ctx));
141:   ctx->st = st;
142:   ctx->alpha = alpha;
143:   ctx->nmat = matIdx?nmat:st->nmat;
144:   PetscCall(PetscMalloc1(ctx->nmat,&ctx->matIdx));
145:   if (matIdx) {
146:     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
147:   } else {
148:     ctx->matIdx[0] = 0;
149:     if (ctx->nmat>1) ctx->matIdx[1] = 1;
150:   }
151:   if (coeffs) {
152:     PetscCall(PetscMalloc1(ctx->nmat,&ctx->coeffs));
153:     for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
154:   }
155:   PetscCall(MatCreateVecs(st->A[0],&ctx->z,NULL));
156:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat));
157:   PetscCall(MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell));
158:   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
159:   PetscCall(MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell));

161:   PetscCall(MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA));
162:   if (st->nmat>1) {
163:     has = hasA;
164:     for (i=1;i<ctx->nmat;i++) {
165:       PetscCall(MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB));
166:       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
167:     }
168:   }
169:   if ((hasA && st->nmat==1) || has) PetscCall(MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }