Actual source code: stshellmat.c

slepc-3.18.3 2023-03-24
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:   MatShellGetContext(A,&ctx);
 31:   ctx->alpha = alpha;
 32:   PetscObjectStateIncrease((PetscObject)A);
 33:   return 0;
 34: }

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

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

 63: static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
 64: {
 65:   ST_MATSHELL    *ctx;
 66:   ST             st;
 67:   PetscInt       i;
 68:   PetscScalar    t=1.0,c;

 70:   MatShellGetContext(A,&ctx);
 71:   st = ctx->st;
 72:   MatMultTranspose(st->A[ctx->matIdx[0]],x,y);
 73:   if (ctx->coeffs && ctx->coeffs[0]!=1.0) VecScale(y,ctx->coeffs[0]);
 74:   if (ctx->alpha!=0.0) {
 75:     for (i=1;i<ctx->nmat;i++) {
 76:       MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z);
 77:       t *= ctx->alpha;
 78:       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
 79:       VecAXPY(y,c,ctx->z);
 80:     }
 81:     if (ctx->nmat==1) VecAXPY(y,ctx->alpha,x); /* y = (A + alpha*I) x */
 82:   }
 83:   return 0;
 84: }

 86: static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
 87: {
 88:   ST_MATSHELL    *ctx;
 89:   ST             st;
 90:   Vec            diagb;
 91:   PetscInt       i;
 92:   PetscScalar    t=1.0,c;

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

114: static PetscErrorCode MatDestroy_Shell(Mat A)
115: {
116:   ST_MATSHELL    *ctx;

118:   MatShellGetContext(A,&ctx);
119:   VecDestroy(&ctx->z);
120:   PetscFree(ctx->matIdx);
121:   PetscFree(ctx->coeffs);
122:   PetscFree(ctx);
123:   return 0;
124: }

126: PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
127: {
128:   PetscInt       n,m,N,M,i;
129:   PetscBool      has=PETSC_FALSE,hasA,hasB;
130:   ST_MATSHELL    *ctx;

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

155:   MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA);
156:   if (st->nmat>1) {
157:     has = hasA;
158:     for (i=1;i<ctx->nmat;i++) {
159:       MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB);
160:       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
161:     }
162:   }
163:   if ((hasA && st->nmat==1) || has) MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);
164:   return 0;
165: }