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          19 : PetscErrorCode STMatShellShift(Mat A,PetscScalar alpha)
      27             : {
      28          19 :   ST_MATSHELL    *ctx;
      29             : 
      30          19 :   PetscFunctionBegin;
      31          19 :   PetscCall(MatShellGetContext(A,&ctx));
      32          19 :   ctx->alpha = alpha;
      33          19 :   PetscCall(PetscObjectStateIncrease((PetscObject)A));
      34          19 :   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       57844 : static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
      42             : {
      43       57844 :   ST_MATSHELL    *ctx;
      44       57844 :   ST             st;
      45       57844 :   PetscInt       i;
      46       57844 :   PetscScalar    t=1.0,c;
      47             : 
      48       57844 :   PetscFunctionBegin;
      49       57844 :   PetscCall(MatShellGetContext(A,&ctx));
      50       57844 :   st = ctx->st;
      51       57844 :   PetscCall(MatMult(st->A[ctx->matIdx[0]],x,y));
      52       57844 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
      53       57844 :   if (ctx->alpha!=0.0) {
      54        1092 :     for (i=1;i<ctx->nmat;i++) {
      55         626 :       PetscCall(MatMult(st->A[ctx->matIdx[i]],x,ctx->z));
      56         626 :       t *= ctx->alpha;
      57         626 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      58         626 :       PetscCall(VecAXPY(y,c,ctx->z));
      59             :     }
      60         466 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
      61             :   }
      62       57844 :   PetscFunctionReturn(PETSC_SUCCESS);
      63             : }
      64             : 
      65         125 : static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
      66             : {
      67         125 :   ST_MATSHELL    *ctx;
      68         125 :   ST             st;
      69         125 :   PetscInt       i;
      70         125 :   PetscScalar    t=1.0,c;
      71             : 
      72         125 :   PetscFunctionBegin;
      73         125 :   PetscCall(MatShellGetContext(A,&ctx));
      74         125 :   st = ctx->st;
      75         125 :   PetscCall(MatMultTranspose(st->A[ctx->matIdx[0]],x,y));
      76         125 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
      77         125 :   if (ctx->alpha!=0.0) {
      78         124 :     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         124 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
      85             :   }
      86         125 :   PetscFunctionReturn(PETSC_SUCCESS);
      87             : }
      88             : 
      89             : #if defined(PETSC_USE_COMPLEX)
      90             : static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A,Vec x,Vec y)
      91             : {
      92             :   ST_MATSHELL    *ctx;
      93             :   ST             st;
      94             :   PetscInt       i;
      95             :   PetscScalar    t=1.0,c;
      96             : 
      97             :   PetscFunctionBegin;
      98             :   PetscCall(MatShellGetContext(A,&ctx));
      99             :   st = ctx->st;
     100             :   PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[0]],x,y));
     101             :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,PetscConj(ctx->coeffs[0])));
     102             :   if (ctx->alpha!=0.0) {
     103             :     for (i=1;i<ctx->nmat;i++) {
     104             :       PetscCall(MatMultHermitianTranspose(st->A[ctx->matIdx[i]],x,ctx->z));
     105             :       t *= ctx->alpha;
     106             :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
     107             :       PetscCall(VecAXPY(y,PetscConj(c),ctx->z));
     108             :     }
     109             :     if (ctx->nmat==1) PetscCall(VecAXPY(y,PetscConj(ctx->alpha),x)); /* y = (A + alpha*I) x */
     110             :   }
     111             :   PetscFunctionReturn(PETSC_SUCCESS);
     112             : }
     113             : #endif
     114             : 
     115          21 : static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
     116             : {
     117          21 :   ST_MATSHELL    *ctx;
     118          21 :   ST             st;
     119          21 :   Vec            diagb;
     120          21 :   PetscInt       i;
     121          21 :   PetscScalar    t=1.0,c;
     122             : 
     123          21 :   PetscFunctionBegin;
     124          21 :   PetscCall(MatShellGetContext(A,&ctx));
     125          21 :   st = ctx->st;
     126          21 :   PetscCall(MatGetDiagonal(st->A[ctx->matIdx[0]],diag));
     127          21 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(diag,ctx->coeffs[0]));
     128          21 :   if (ctx->alpha!=0.0) {
     129          18 :     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          21 :   PetscFunctionReturn(PETSC_SUCCESS);
     142             : }
     143             : 
     144          82 : static PetscErrorCode MatDestroy_Shell(Mat A)
     145             : {
     146          82 :   ST_MATSHELL    *ctx;
     147             : 
     148          82 :   PetscFunctionBegin;
     149          82 :   PetscCall(MatShellGetContext(A,&ctx));
     150          82 :   PetscCall(VecDestroy(&ctx->z));
     151          82 :   PetscCall(PetscFree(ctx->matIdx));
     152          82 :   PetscCall(PetscFree(ctx->coeffs));
     153          82 :   PetscCall(PetscFree(ctx));
     154          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     155             : }
     156             : 
     157          82 : PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
     158             : {
     159          82 :   PetscInt       n,m,N,M,i;
     160          82 :   PetscBool      has=PETSC_FALSE,hasA,hasB;
     161          82 :   ST_MATSHELL    *ctx;
     162             : 
     163          82 :   PetscFunctionBegin;
     164          82 :   PetscCall(MatGetSize(st->A[0],&M,&N));
     165          82 :   PetscCall(MatGetLocalSize(st->A[0],&m,&n));
     166          82 :   PetscCall(PetscNew(&ctx));
     167          82 :   ctx->st = st;
     168          82 :   ctx->alpha = alpha;
     169          82 :   ctx->nmat = matIdx?nmat:st->nmat;
     170          82 :   PetscCall(PetscMalloc1(ctx->nmat,&ctx->matIdx));
     171          82 :   if (matIdx) {
     172          39 :     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
     173             :   } else {
     174          72 :     ctx->matIdx[0] = 0;
     175          72 :     if (ctx->nmat>1) ctx->matIdx[1] = 1;
     176             :   }
     177          82 :   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          82 :   PetscCall(MatCreateVecs(st->A[0],&ctx->z,NULL));
     182          82 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat));
     183          82 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell));
     184          82 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
     185             : #if defined(PETSC_USE_COMPLEX)
     186             :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_Shell));
     187             : #else
     188          82 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
     189             : #endif
     190          82 :   PetscCall(MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell));
     191             : 
     192          82 :   PetscCall(MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA));
     193          82 :   if (st->nmat>1) {
     194          35 :     has = hasA;
     195          79 :     for (i=1;i<ctx->nmat;i++) {
     196          44 :       PetscCall(MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB));
     197          44 :       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
     198             :     }
     199             :   }
     200          82 :   if ((hasA && st->nmat==1) || has) PetscCall(MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
     201          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     202             : }