LCOV - code coverage report
Current view: top level - sys/classes/st/interface - stshellmat.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 92 98 93.9 %
Date: 2020-05-28 03:14:05 Functions: 6 6 100.0 %

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

Generated by: LCOV version 1.13