LCOV - code coverage report
Current view: top level - sys/classes/st/interface - stshellmat.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 91 97 93.8 %
Date: 2019-09-16 06:01:52 Functions: 6 6 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, 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_SHELLMAT;
      25             : 
      26          18 : PetscErrorCode STMatShellShift(Mat A,PetscScalar alpha)
      27             : {
      28             :   PetscErrorCode ierr;
      29             :   ST_SHELLMAT    *ctx;
      30             : 
      31          18 :   PetscFunctionBegin;
      32          18 :   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
      33          18 :   ctx->alpha = alpha;
      34          18 :   PetscFunctionReturn(0);
      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        1784 : static PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
      42             : {
      43             :   PetscErrorCode ierr;
      44             :   ST_SHELLMAT    *ctx;
      45             :   ST             st;
      46             :   PetscInt       i;
      47        1784 :   PetscScalar    t=1.0,c;
      48             : 
      49        1784 :   PetscFunctionBegin;
      50        1784 :   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
      51        1784 :   st = ctx->st;
      52        1784 :   ierr = MatMult(st->A[ctx->matIdx[0]],x,y);CHKERRQ(ierr);
      53        1784 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
      54         182 :     ierr = VecScale(y,ctx->coeffs[0]);CHKERRQ(ierr);
      55             :   }
      56        1784 :   if (ctx->alpha!=0.0) {
      57         671 :     for (i=1;i<ctx->nmat;i++) {
      58         671 :       ierr = MatMult(st->A[ctx->matIdx[i]],x,ctx->z);CHKERRQ(ierr);
      59         671 :       t *= ctx->alpha;
      60         671 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      61         671 :       ierr = VecAXPY(y,c,ctx->z);CHKERRQ(ierr);
      62             :     }
      63         483 :     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
      64         159 :       ierr = VecAXPY(y,ctx->alpha,x);CHKERRQ(ierr);
      65             :     }
      66             :   }
      67        1784 :   PetscFunctionReturn(0);
      68             : }
      69             : 
      70         119 : static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
      71             : {
      72             :   PetscErrorCode ierr;
      73             :   ST_SHELLMAT    *ctx;
      74             :   ST             st;
      75             :   PetscInt       i;
      76         119 :   PetscScalar    t=1.0,c;
      77             : 
      78         119 :   PetscFunctionBegin;
      79         119 :   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
      80         119 :   st = ctx->st;
      81         119 :   ierr = MatMultTranspose(st->A[ctx->matIdx[0]],x,y);CHKERRQ(ierr);
      82         119 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
      83           0 :     ierr = VecScale(y,ctx->coeffs[0]);CHKERRQ(ierr);
      84             :   }
      85         119 :   if (ctx->alpha!=0.0) {
      86           0 :     for (i=1;i<ctx->nmat;i++) {
      87           0 :       ierr = MatMultTranspose(st->A[ctx->matIdx[i]],x,ctx->z);CHKERRQ(ierr);
      88           0 :       t *= ctx->alpha;
      89           0 :       c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
      90           0 :       ierr = VecAXPY(y,c,ctx->z);CHKERRQ(ierr);
      91             :     }
      92         118 :     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
      93         118 :       ierr = VecAXPY(y,ctx->alpha,x);CHKERRQ(ierr);
      94             :     }
      95             :   }
      96         119 :   PetscFunctionReturn(0);
      97             : }
      98             : 
      99          10 : static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec diag)
     100             : {
     101             :   PetscErrorCode ierr;
     102             :   ST_SHELLMAT    *ctx;
     103             :   ST             st;
     104             :   Vec            diagb;
     105             :   PetscInt       i;
     106          10 :   PetscScalar    t=1.0,c;
     107             : 
     108          10 :   PetscFunctionBegin;
     109          10 :   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
     110          10 :   st = ctx->st;
     111          10 :   ierr = MatGetDiagonal(st->A[ctx->matIdx[0]],diag);CHKERRQ(ierr);
     112          10 :   if (ctx->coeffs && ctx->coeffs[0]!=1.0) {
     113           2 :     ierr = VecScale(diag,ctx->coeffs[0]);CHKERRQ(ierr);
     114             :   }
     115          10 :   if (ctx->alpha!=0.0) {
     116          10 :     if (ctx->nmat==1) {    /* y = (A + alpha*I) x */
     117           5 :       ierr = VecShift(diag,ctx->alpha);CHKERRQ(ierr);
     118             :     } else {
     119           5 :       ierr = VecDuplicate(diag,&diagb);CHKERRQ(ierr);
     120           9 :       for (i=1;i<ctx->nmat;i++) {
     121           9 :         ierr = MatGetDiagonal(st->A[ctx->matIdx[i]],diagb);CHKERRQ(ierr);
     122           9 :         t *= ctx->alpha;
     123           9 :         c = (ctx->coeffs)?t*ctx->coeffs[i]:t;
     124           9 :         ierr = VecAYPX(diag,c,diagb);CHKERRQ(ierr);
     125             :       }
     126           5 :       ierr = VecDestroy(&diagb);CHKERRQ(ierr);
     127             :     }
     128             :   }
     129          10 :   PetscFunctionReturn(0);
     130             : }
     131             : 
     132          72 : static PetscErrorCode MatDestroy_Shell(Mat A)
     133             : {
     134             :   PetscErrorCode ierr;
     135             :   ST_SHELLMAT    *ctx;
     136             : 
     137          72 :   PetscFunctionBegin;
     138          72 :   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
     139          72 :   ierr = VecDestroy(&ctx->z);CHKERRQ(ierr);
     140          72 :   ierr = PetscFree(ctx->matIdx);CHKERRQ(ierr);
     141          72 :   ierr = PetscFree(ctx->coeffs);CHKERRQ(ierr);
     142          72 :   ierr = PetscFree(ctx);CHKERRQ(ierr);
     143          72 :   PetscFunctionReturn(0);
     144             : }
     145             : 
     146          72 : PetscErrorCode STMatShellCreate(ST st,PetscScalar alpha,PetscInt nmat,PetscInt *matIdx,PetscScalar *coeffs,Mat *mat)
     147             : {
     148             :   PetscErrorCode ierr;
     149             :   PetscInt       n,m,N,M,i;
     150          72 :   PetscBool      has=PETSC_FALSE,hasA,hasB;
     151             :   ST_SHELLMAT    *ctx;
     152             : 
     153          72 :   PetscFunctionBegin;
     154          72 :   ierr = MatGetSize(st->A[0],&M,&N);CHKERRQ(ierr);
     155          72 :   ierr = MatGetLocalSize(st->A[0],&m,&n);CHKERRQ(ierr);
     156          72 :   ierr = PetscNew(&ctx);CHKERRQ(ierr);
     157          72 :   ctx->st = st;
     158          72 :   ctx->alpha = alpha;
     159          72 :   ctx->nmat = matIdx?nmat:st->nmat;
     160          72 :   ierr = PetscMalloc1(ctx->nmat,&ctx->matIdx);CHKERRQ(ierr);
     161          72 :   if (matIdx) {
     162          29 :     for (i=0;i<ctx->nmat;i++) ctx->matIdx[i] = matIdx[i];
     163             :   } else {
     164          62 :     ctx->matIdx[0] = 0;
     165          62 :     if (ctx->nmat>1) ctx->matIdx[1] = 1;
     166             :   }
     167          72 :   if (coeffs) {
     168          10 :     ierr = PetscMalloc1(ctx->nmat,&ctx->coeffs);CHKERRQ(ierr);
     169          29 :     for (i=0;i<ctx->nmat;i++) ctx->coeffs[i] = coeffs[i];
     170             :   }
     171          72 :   ierr = MatCreateVecs(st->A[0],&ctx->z,NULL);CHKERRQ(ierr);
     172          72 :   ierr = MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,(void*)ctx,mat);CHKERRQ(ierr);
     173          72 :   ierr = MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatMult_Shell);CHKERRQ(ierr);
     174          72 :   ierr = MatShellSetOperation(*mat,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell);CHKERRQ(ierr);
     175          72 :   ierr = MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell);CHKERRQ(ierr);
     176             : 
     177          72 :   ierr = MatHasOperation(st->A[ctx->matIdx[0]],MATOP_GET_DIAGONAL,&hasA);CHKERRQ(ierr);
     178          72 :   if (st->nmat>1) {
     179          14 :     has = hasA;
     180          37 :     for (i=1;i<ctx->nmat;i++) {
     181          23 :       ierr = MatHasOperation(st->A[ctx->matIdx[i]],MATOP_GET_DIAGONAL,&hasB);CHKERRQ(ierr);
     182          23 :       has = (has && hasB)? PETSC_TRUE: PETSC_FALSE;
     183             :     }
     184             :   }
     185          72 :   if ((hasA && st->nmat==1) || has) {
     186          71 :     ierr = MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);CHKERRQ(ierr);
     187             :   }
     188          72 :   PetscFunctionReturn(0);
     189             : }
     190             : 

Generated by: LCOV version 1.13