LCOV - code coverage report
Current view: top level - sys/classes/st/interface - stshellmat.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 116 124 93.5 %
Date: 2024-03-28 00:28:38 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

          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        4759 : static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
      42             : {
      43        4759 :   ST_MATSHELL    *ctx;
      44        4759 :   ST             st;
      45        4759 :   PetscInt       i;
      46        4759 :   PetscScalar    t=1.0,c;
      47             : 
      48        4759 :   PetscFunctionBegin;
      49        4759 :   PetscCall(MatShellGetContext(A,&ctx));
      50        4759 :   st = ctx->st;
      51        4759 :   PetscCall(MatMult(st->A[ctx->matIdx[0]],x,y));
      52        4759 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) PetscCall(VecScale(y,ctx->coeffs[0]));
      53        4759 :   if (ctx->alpha!=0.0) {
      54        1305 :     for (i=1;i<ctx->nmat;i++) {
      55         648 :       PetscCall(MatMult(st->A[ctx->matIdx[i]],x,ctx->z));
      56         648 :       t *= ctx->alpha;
      57         648 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      58         648 :       PetscCall(VecAXPY(y,c,ctx->z));
      59             :     }
      60         657 :     if (ctx->nmat==1) PetscCall(VecAXPY(y,ctx->alpha,x)); /* y = (A + alpha*I) x */
      61             :   }
      62        4759 :   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          85 : static PetscErrorCode MatDestroy_Shell(Mat A)
     145             : {
     146          85 :   ST_MATSHELL    *ctx;
     147             : 
     148          85 :   PetscFunctionBegin;
     149          85 :   PetscCall(MatShellGetContext(A,&ctx));
     150          85 :   PetscCall(VecDestroy(&ctx->z));
     151          85 :   PetscCall(PetscFree(ctx->matIdx));
     152          85 :   PetscCall(PetscFree(ctx->coeffs));
     153          85 :   PetscCall(PetscFree(ctx));
     154          85 :   PetscFunctionReturn(PETSC_SUCCESS);
     155             : }
     156             : 
     157          85 : PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
     158             : {
     159          85 :   PetscInt       n,m,N,M,i;
     160          85 :   PetscBool      has=PETSC_FALSE,hasA,hasB;
     161          85 :   ST_MATSHELL    *ctx;
     162             : 
     163          85 :   PetscFunctionBegin;
     164          85 :   PetscCall(MatGetSize(st->A[0],&M,&N));
     165          85 :   PetscCall(MatGetLocalSize(st->A[0],&m,&n));
     166          85 :   PetscCall(PetscNew(&ctx));
     167          85 :   ctx->st = st;
     168          85 :   ctx->alpha = alpha;
     169          85 :   ctx->nmat = matIdx?nmat:st->nmat;
     170          85 :   PetscCall(PetscMalloc1(ctx->nmat,&ctx->matIdx));
     171          85 :   if (matIdx) {
     172          39 :     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
     173             :   } else {
     174          75 :     ctx->matIdx[0] = 0;
     175          75 :     if (ctx->nmat>1) ctx->matIdx[1] = 1;
     176             :   }
     177          85 :   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          85 :   PetscCall(MatCreateVecs(st->A[0],&ctx->z,NULL));
     182          85 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat));
     183          85 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell));
     184          85 :   PetscCall(MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
     185             : #if defined(PETSC_USE_COMPLEX)
     186          85 :   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          85 :   PetscCall(MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell));
     191             : 
     192          85 :   PetscCall(MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA));
     193          85 :   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          85 :   if ((hasA && st->nmat==1) || has) PetscCall(MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
     201          85 :   PetscFunctionReturn(PETSC_SUCCESS);
     202             : }

Generated by: LCOV version 1.14