Actual source code: stshellmat.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }