Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      14             : 
      15             : #include <slepc/private/stimpl.h>
      16             : 
      17             : typedef struct {
      18             :   PetscScalar alpha;
      19             :   PetscScalar *coeffs;
      20             :   ST          st;
      21             :   Vec         z;
      22             :   PetscInt    nmat;
      23             :   PetscInt    *matIdx;
      24             : } ST_MATSHELL;
      25             : 
      26          24 : PetscErrorCode STMatShellShift(Mat A,PetscScalar alpha)
      27             : {
      28          24 :   ST_MATSHELL    *ctx;
      29             : 
      30          24 :   PetscFunctionBegin;
      31          24 :   PetscCall(MatShellGetContext(A,&ctx));
      32          24 :   ctx->alpha = alpha;
      33          24 :   PetscCall(PetscObjectStateIncrease((PetscObject)A));
      34          24 :   PetscFunctionReturn(PETSC_SUCCESS);
      35             : }
      36             : 
      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        4717 : static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
      42             : {
      43        4717 :   ST_MATSHELL    *ctx;
      44        4717 :   ST             st;
      45        4717 :   PetscInt       i;
      46        4717 :   PetscScalar    t=1.0,c;
      47             : 
      48        4717 :   PetscFunctionBegin;
      49        4717 :   PetscCall(MatShellGetContext(A,&ctx));
      50        4717 :   st = ctx->st;
      51        4717 :   PetscCall(MatMult(st->A[ctx->matIdx[0]],x,y));
      52        4717 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
      53        4717 :   if (ctx->alpha!=0.0) {
      54        1301 :     for (i=1;i<ctx->nmat;i++) {
      55         645 :       PetscCall(MatMult(st->A[ctx->matIdx[i]],x,ctx->z));
      56         645 :       t *= ctx->alpha;
      57         645 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      58         645 :       PetscCall(VecAXPY(y,c,ctx->z));
      59             :     }
      60         656 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
      61             :   }
      62        4717 :   PetscFunctionReturn(PETSC_SUCCESS);
      63             : }
      64             : 
      65         354 : static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
      66             : {
      67         354 :   ST_MATSHELL    *ctx;
      68         354 :   ST             st;
      69         354 :   PetscInt       i;
      70         354 :   PetscScalar    t=1.0,c;
      71             : 
      72         354 :   PetscFunctionBegin;
      73         354 :   PetscCall(MatShellGetContext(A,&ctx));
      74         354 :   st = ctx->st;
      75         354 :   PetscCall(MatMultTranspose(st->A[ctx->matIdx[0]],x,y));
      76         354 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
      77         354 :   if (ctx->alpha!=0.0) {
      78         354 :     for (i=1;i<ctx->nmat;i++) {
      79           0 :       PetscCall(MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z));
      80           0 :       t *= ctx->alpha;
      81           0 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      82           0 :       PetscCall(VecAXPY(y,c,ctx->z));
      83             :     }
      84         354 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
      85             :   }
      86         354 :   PetscFunctionReturn(PETSC_SUCCESS);
      87             : }
      88             : 
      89             : #if defined(PETSC_USE_COMPLEX)
      90           5 : static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A,Vec x,Vec y)
      91             : {
      92           5 :   ST_MATSHELL    *ctx;
      93           5 :   ST             st;
      94           5 :   PetscInt       i;
      95           5 :   PetscScalar    t=1.0,c;
      96             : 
      97           5 :   PetscFunctionBegin;
      98           5 :   PetscCall(MatShellGetContext(A,&ctx));
      99           5 :   st = ctx->st;
     100           5 :   PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[0]],x,y));
     101           5 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,PetscConj(ctx->coeffs[0])));
     102           5 :   if (ctx->alpha!=0.0) {
     103           4 :     for (i=1;i<ctx->nmat;i++) {
     104           0 :       PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[i]],x,ctx->z));
     105           0 :       t *= ctx->alpha;
     106           0 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
     107           0 :       PetscCall(VecAXPY(y,PetscConj(c),ctx->z));
     108             :     }
     109           4 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,PetscConj(ctx->alpha),x)); /* y = (A + alpha*I) x */
     110             :   }
     111           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     112             : }
     113             : #endif
     114             : 
     115          26 : static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
     116             : {
     117          26 :   ST_MATSHELL    *ctx;
     118          26 :   ST             st;
     119          26 :   Vec            diagb;
     120          26 :   PetscInt       i;
     121          26 :   PetscScalar    t=1.0,c;
     122             : 
     123          26 :   PetscFunctionBegin;
     124          26 :   PetscCall(MatShellGetContext(A,&ctx));
     125          26 :   st = ctx->st;
     126          26 :   PetscCall(MatGetDiagonal(st->A[ctx->matIdx[0]],diag));
     127          26 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(diag,ctx->coeffs[0]));
     128          26 :   if (ctx->alpha!=0.0) {
     129          24 :     if (ctx->nmat==1) PetscCall(VecShift(diag,ctx->alpha)); /* y = (A + alpha*I) x */
     130             :     else {
     131           8 :       PetscCall(VecDuplicate(diag,&diagb));
     132          22 :       for (i=1;i<ctx->nmat;i++) {
     133          14 :         PetscCall(MatGetDiagonal(st->A[ctx->matIdx[i]],diagb));
     134          14 :         t *= ctx->alpha;
     135          14 :         c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
     136          14 :         PetscCall(VecAYPX(diag,c,diagb));
     137             :       }
     138           8 :       PetscCall(VecDestroy(&diagb));
     139             :     }
     140             :   }
     141          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     142             : }
     143             : 
     144          83 : static PetscErrorCode MatDestroy_Shell(Mat A)
     145             : {
     146          83 :   ST_MATSHELL    *ctx;
     147             : 
     148          83 :   PetscFunctionBegin;
     149          83 :   PetscCall(MatShellGetContext(A,&ctx));
     150          83 :   PetscCall(VecDestroy(&ctx->z));
     151          83 :   PetscCall(PetscFree(ctx->matIdx));
     152          83 :   PetscCall(PetscFree(ctx->coeffs));
     153          83 :   PetscCall(PetscFree(ctx));
     154          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     155             : }
     156             : 
     157          83 : PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
     158             : {
     159          83 :   PetscInt       n,m,N,M,i;
     160          83 :   PetscBool      has=PETSC_FALSE,hasA,hasB;
     161          83 :   ST_MATSHELL    *ctx;
     162             : 
     163          83 :   PetscFunctionBegin;
     164          83 :   PetscCall(MatGetSize(st->A[0],&M,&N));
     165          83 :   PetscCall(MatGetLocalSize(st->A[0],&m,&n));
     166          83 :   PetscCall(PetscNew(&ctx));
     167          83 :   ctx->st = st;
     168          83 :   ctx->alpha = alpha;
     169          83 :   ctx->nmat = matIdx?nmat:st->nmat;
     170          83 :   PetscCall(PetscMalloc1(ctx->nmat,&ctx->matIdx));
     171          83 :   if (matIdx) {
     172          39 :     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
     173             :   } else {
     174          73 :     ctx->matIdx[0] = 0;
     175          73 :     if (ctx->nmat>1) ctx->matIdx[1] = 1;
     176             :   }
     177          83 :   if (coeffs) {
     178          10 :     PetscCall(PetscMalloc1(ctx->nmat,&ctx->coeffs));
     179          39 :     for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
     180             :   }
     181          83 :   PetscCall(MatCreateVecs(st->A[0],&ctx->z,NULL));
     182          83 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat));
     183          83 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell));
     184          83 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
     185             : #if defined(PETSC_USE_COMPLEX)
     186          83 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_Shell));
     187             : #else
     188             :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
     189             : #endif
     190          83 :   PetscCall(MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell));
     191             : 
     192          83 :   PetscCall(MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA));
     193          83 :   if (st->nmat>1) {
     194          31 :     has = hasA;
     195          71 :     for (i=1;i<ctx->nmat;i++) {
     196          40 :       PetscCall(MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB));
     197          40 :       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
     198             :     }
     199             :   }
     200          83 :   if ((hasA && st->nmat==1) || has) PetscCall(MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
     201          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     202             : }