Actual source code: stshellmat.c

slepc-3.13.2 2020-05-12
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: {
 29:   ST_MATSHELL    *ctx;

 32:   MatShellGetContext(A,(void**)&ctx);
 33:   ctx->alpha = alpha;
 34:   PetscObjectStateIncrease((PetscObject)A);
 35:   return(0);
 36: }

 38: /*
 39:   For i=0:nmat-1 computes y = (sum_i (coeffs[i]*alpha^i*st->A[idx[i]]))x
 40:   If null coeffs computes with coeffs[i]=1.0
 41: */
 42: static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
 43: {
 45:   ST_MATSHELL    *ctx;
 46:   ST             st;
 47:   PetscInt       i;
 48:   PetscScalar    t=1.0,c;

 51:   MatShellGetContext(A,(void**)&ctx);
 52:   st = ctx->st;
 53:   MatMult(st->A[ctx->matIdx[0]],x,y);
 54:   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
 55:     VecScale(y,ctx->coeffs[0]);
 56:   }
 57:   if (ctx->alpha!=0.0) {
 58:     for (i=1;i<ctx->nmat;i++) {
 59:       MatMult(st->A[ctx->matIdx[i]],x,ctx->z);
 60:       t *= ctx->alpha;
 61:       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
 62:       VecAXPY(y,c,ctx->z);
 63:     }
 64:     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
 65:       VecAXPY(y,ctx->alpha,x);
 66:     }
 67:   }
 68:   return(0);
 69: }

 71: static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
 72: {
 74:   ST_MATSHELL    *ctx;
 75:   ST             st;
 76:   PetscInt       i;
 77:   PetscScalar    t=1.0,c;

 80:   MatShellGetContext(A,(void**)&ctx);
 81:   st = ctx->st;
 82:   MatMultTranspose(st->A[ctx->matIdx[0]],x,y);
 83:   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
 84:     VecScale(y,ctx->coeffs[0]);
 85:   }
 86:   if (ctx->alpha!=0.0) {
 87:     for (i=1;i<ctx->nmat;i++) {
 88:       MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z);
 89:       t *= ctx->alpha;
 90:       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
 91:       VecAXPY(y,c,ctx->z);
 92:     }
 93:     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
 94:       VecAXPY(y,ctx->alpha,x);
 95:     }
 96:   }
 97:   return(0);
 98: }

100: static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
101: {
103:   ST_MATSHELL    *ctx;
104:   ST             st;
105:   Vec            diagb;
106:   PetscInt       i;
107:   PetscScalar    t=1.0,c;

110:   MatShellGetContext(A,(void**)&ctx);
111:   st = ctx->st;
112:   MatGetDiagonal(st->A[ctx->matIdx[0]],diag);
113:   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
114:     VecScale(diag,ctx->coeffs[0]);
115:   }
116:   if (ctx->alpha!=0.0) {
117:     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
118:       VecShift(diag,ctx->alpha);
119:     } else {
120:       VecDuplicate(diag,&diagb);
121:       for (i=1;i<ctx->nmat;i++) {
122:         MatGetDiagonal(st->A[ctx->matIdx[i]],diagb);
123:         t *= ctx->alpha;
124:         c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
125:         VecAYPX(diag,c,diagb);
126:       }
127:       VecDestroy(&diagb);
128:     }
129:   }
130:   return(0);
131: }

133: static PetscErrorCode MatDestroy_Shell(Mat A)
134: {
136:   ST_MATSHELL    *ctx;

139:   MatShellGetContext(A,(void**)&ctx);
140:   VecDestroy(&ctx->z);
141:   PetscFree(ctx->matIdx);
142:   PetscFree(ctx->coeffs);
143:   PetscFree(ctx);
144:   return(0);
145: }

147: PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
148: {
150:   PetscInt       n,m,N,M,i;
151:   PetscBool      has=PETSC_FALSE,hasA,hasB;
152:   ST_MATSHELL    *ctx;

155:   MatGetSize(st->A[0],&M,&N);
156:   MatGetLocalSize(st->A[0],&m,&n);
157:   PetscNew(&ctx);
158:   ctx->st = st;
159:   ctx->alpha = alpha;
160:   ctx->nmat = matIdx?nmat:st->nmat;
161:   PetscMalloc1(ctx->nmat,&ctx->matIdx);
162:   if (matIdx) {
163:     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
164:   } else {
165:     ctx->matIdx[0] = 0;
166:     if (ctx->nmat>1) ctx->matIdx[1] = 1;
167:   }
168:   if (coeffs) {
169:     PetscMalloc1(ctx->nmat,&ctx->coeffs);
170:     for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
171:   }
172:   MatCreateVecs(st->A[0],&ctx->z,NULL);
173:   MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat);
174:   MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell);
175:   MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell);
176:   MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell);

178:   MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA);
179:   if (st->nmat>1) {
180:     has = hasA;
181:     for (i=1;i<ctx->nmat;i++) {
182:       MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB);
183:       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
184:     }
185:   }
186:   if ((hasA && st->nmat==1) || has) {
187:     MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);
188:   }
189:   return(0);
190: }