LCOV - code coverage report
Current view: top level - lme/interface - lmedense.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 89 89 100.0 %
Date: 2019-09-16 06:01:52 Functions: 5 5 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             :    Routines for solving dense matrix equations, in some cases calling SLICOT
      12             : */
      13             : 
      14             : #include <slepc/private/lmeimpl.h>     /*I "slepclme.h" I*/
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : /*
      18             :    LMERankSVD - given a square matrix L, compute its SVD U*S*V', and determine the
      19             :    numerical rank. On exit, U contains U*S and L is overwritten with V'
      20             : */
      21          10 : PetscErrorCode LMERankSVD(LME lme,PetscInt n,PetscScalar *L,PetscScalar *U,PetscInt *rank)
      22             : {
      23             : #if defined(PETSC_MISSING_LAPACK_GESVD)
      24             :   PetscFunctionBegin;
      25             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable");
      26             : #else
      27             :   PetscErrorCode ierr;
      28          10 :   PetscInt       i,j,rk=0;
      29             :   PetscScalar    *work;
      30             :   PetscReal      tol,*sg,*rwork;
      31             :   PetscBLASInt   n_,info,lw_;
      32             : 
      33          10 :   PetscFunctionBegin;
      34          10 :   ierr = PetscCalloc3(n,&sg,10*n,&work,5*n,&rwork);CHKERRQ(ierr);
      35          10 :   ierr = PetscBLASIntCast(n,&n_);CHKERRQ(ierr);
      36          10 :   lw_ = 10*n_;
      37             : #if !defined (PETSC_USE_COMPLEX)
      38          10 :   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,L,&n_,sg,U,&n_,NULL,&n_,work,&lw_,&info));
      39             : #else
      40             :   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,L,&n_,sg,U,&n_,NULL,&n_,work,&lw_,rwork,&info));
      41             : #endif
      42          10 :   SlepcCheckLapackInfo("gesvd",info);
      43          10 :   tol = 10*PETSC_MACHINE_EPSILON*n*sg[0];
      44         382 :   for (j=0;j<n;j++) {
      45         372 :     if (sg[j]>tol) {
      46       14184 :       for (i=0;i<n;i++) U[i+j*n] *= sg[j];
      47         372 :       rk++;
      48             :     } else break;
      49             :   }
      50          10 :   *rank = rk;
      51          10 :   ierr = PetscFree3(sg,work,rwork);CHKERRQ(ierr);
      52          10 :   PetscFunctionReturn(0);
      53             : #endif
      54             : }
      55             : 
      56             : #if defined(PETSC_USE_INFO)
      57             : /*
      58             :    LyapunovResidual - compute the residual norm ||H*L*L'+L*L'*H'+r*r'||
      59             : */
      60           4 : static PetscErrorCode LyapunovResidual(PetscScalar *H,PetscInt m,PetscInt ldh,PetscScalar *r,PetscScalar *L,PetscInt ldl,PetscReal *res)
      61             : {
      62             :   PetscErrorCode ierr;
      63             :   PetscBLASInt   n,ld;
      64             :   PetscInt       i,j;
      65           4 :   PetscScalar    *M,*R,zero=0.0,done=1.0;
      66             : 
      67           4 :   PetscFunctionBegin;
      68           4 :   *res = 0;
      69           4 :   ierr = PetscBLASIntCast(ldh,&ld);CHKERRQ(ierr);
      70           4 :   ierr = PetscBLASIntCast(m,&n);CHKERRQ(ierr);
      71           4 :   ierr = PetscMalloc2(m*m,&M,m*m,&R);CHKERRQ(ierr);
      72             : 
      73             :   /* R = r*r' */
      74         126 :   for (i=0;i<m;i++) {
      75        4410 :     for (j=0;j<m;j++) R[i+j*m] = r[i]*r[j];
      76             :   }
      77             :   /* M = H*L */
      78           4 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,H,&ld,L,&n,&zero,M,&n));
      79             :   /* R = R+M*L' */
      80           4 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,M,&n,L,&n,&done,R,&n));
      81             :   /* R = R+L*M' */
      82           4 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,L,&n,M,&n,&done,R,&n));
      83             :   
      84           4 :   *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
      85           4 :   ierr = PetscFree2(M,R);CHKERRQ(ierr);
      86           4 :   PetscFunctionReturn(0);
      87             : }
      88             : #endif
      89             : 
      90             : #if defined(SLEPC_HAVE_SLICOT)
      91             : /*
      92             :    LyapunovFact_SLICOT - implementation used when SLICOT is available
      93             : */
      94             : static PetscErrorCode LyapunovChol_SLICOT(PetscScalar *H,PetscInt m,PetscInt ldh,PetscScalar *r,PetscScalar *L,PetscInt ldl,PetscReal *res)
      95             : {
      96             : #if defined(PETSC_MISSING_LAPACK_HSEQR)
      97             :   PetscFunctionBegin;
      98             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"HSEQR - Lapack routine is unavailable");
      99             : #else
     100             :   PetscErrorCode ierr;
     101             :   PetscBLASInt   ilo=1,lwork,info,n,ld,ld1,ione=1;
     102             :   PetscInt       i,j;
     103             :   PetscReal      scal;
     104             :   PetscScalar    *Q,*W,*z,*wr,*work,zero=0.0,done=1.0,alpha,beta;
     105             : #if !defined(PETSC_USE_COMPLEX)
     106             :   PetscScalar    *wi;
     107             : #endif
     108             : 
     109             :   PetscFunctionBegin;
     110             :   ierr = PetscBLASIntCast(ldh,&ld);CHKERRQ(ierr);
     111             :   ierr = PetscBLASIntCast(ldl,&ld1);CHKERRQ(ierr);
     112             :   ierr = PetscBLASIntCast(m,&n);CHKERRQ(ierr);
     113             :   ierr = PetscBLASIntCast(5*m,&lwork);CHKERRQ(ierr);
     114             : 
     115             :   /* compute the (real) Schur form of H */
     116             : #if !defined(PETSC_USE_COMPLEX)
     117             :   ierr = PetscCalloc6(m*m,&Q,m*m,&W,m,&z,m,&wr,m,&wi,5*m,&work);CHKERRQ(ierr);
     118             :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,H,&ld,wr,wi,Q,&n,work,&lwork,&info));
     119             : #else
     120             :   ierr = PetscCalloc5(m*m,&Q,m*m,&W,m,&z,m,&wr,5*m,&work);CHKERRQ(ierr);
     121             :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,H,&ld,wr,Q,&n,work,&lwork,&info));
     122             : #endif
     123             :   SlepcCheckLapackInfo("hseqr",info);
     124             : #if defined(PETSC_USE_DEBUG)
     125             :   for (i=0;i<m;i++) if (PetscRealPart(wr[i])>0.0) SETERRQ(PETSC_COMM_SELF,1,"Positive eigenvalue found, the coefficient matrix is not stable");
     126             : #endif
     127             : 
     128             :   /* copy r into first column of W */
     129             :   ierr = PetscArraycpy(W,r,m);CHKERRQ(ierr);
     130             : 
     131             :   /* solve Lyapunov equation (Hammarling) */
     132             :   PetscStackCallBLAS("SLICOTsb03od",SLICOTsb03od_("C","F","N",&n,&ione,H,&ld,Q,&n,W,&n,&scal,wr,wi,work,&lwork,&info));
     133             :   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SLICOT subroutine SB03OD %d",(int)info);
     134             :   if (scal!=1.0) SETERRQ1(PETSC_COMM_SELF,1,"Current implementation cannot handle scale factor %g",scal);
     135             : 
     136             :   /* Tranpose L = W' */
     137             :   for (j=0;j<m;j++) {
     138             :     for (i=j;i<m;i++) L[i+j*ldl] = W[j+i*m];
     139             :   }
     140             : 
     141             :   /* resnorm = norm(H(m+1,:)*L*L'), use z = L*L(m,:)' */
     142             :   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&done,L,&ld1,W+(m-1)*m,&ione,&zero,z,&ione));
     143             :   *res = 0.0;
     144             :   beta = H[m+(m-1)*ldh];
     145             :   for (j=0;j<m;j++) {
     146             :     alpha = beta*z[j];
     147             :     *res += alpha*alpha;
     148             :   }
     149             :   *res = PetscSqrtReal(*res);
     150             : 
     151             : #if !defined(PETSC_USE_COMPLEX)
     152             :   ierr = PetscFree6(Q,W,z,wr,wi,work);CHKERRQ(ierr);
     153             : #else
     154             :   ierr = PetscFree5(Q,W,z,wr,work);CHKERRQ(ierr);
     155             : #endif
     156             :   PetscFunctionReturn(0);
     157             : #endif
     158             : }
     159             : 
     160             : #else
     161             : 
     162             : #if 0
     163             : /*
     164             :    AbsEig - given a matrix A that may be slightly indefinite (hence Cholesky fails)
     165             :    modify it by taking the absolute value of the eigenvalues: [U,S] = eig(A); A = U*abs(S)*U';
     166             : */
     167             : static PetscErrorCode AbsEig(PetscScalar *A,PetscInt m)
     168             : {
     169             : #if defined(PETSC_MISSING_LAPACK_SYEV) || defined(SLEPC_MISSING_LAPACK_LACPY)
     170             :   PetscFunctionBegin;
     171             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYEV/LACPY - Lapack routines are unavailable");
     172             : #else
     173             :   PetscErrorCode ierr;
     174             :   PetscInt       i,j;
     175             :   PetscBLASInt   n,ld,lwork,info;
     176             :   PetscScalar    *Q,*W,*work,a,one=1.0,zero=0.0;
     177             :   PetscReal      *eig,dummy;
     178             : #if defined(PETSC_USE_COMPLEX)
     179             :   PetscReal      *rwork,rdummy;
     180             : #endif
     181             : 
     182             :   PetscFunctionBegin;
     183             :   ierr = PetscBLASIntCast(m,&n);CHKERRQ(ierr);
     184             :   ld = n;
     185             : 
     186             :   /* workspace query and memory allocation */
     187             :   lwork = -1;
     188             : #if defined(PETSC_USE_COMPLEX)
     189             :   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,A,&ld,&dummy,&a,&lwork,&rdummy,&info));
     190             :   ierr = PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);CHKERRQ(ierr);
     191             :   ierr = PetscMalloc5(m,&eig,m*m,&Q,m*n,&W,lwork,&work,PetscMax(1,3*m-2),&rwork);CHKERRQ(ierr);
     192             : #else
     193             :   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,A,&ld,&dummy,&a,&lwork,&info));
     194             :   ierr = PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);CHKERRQ(ierr);
     195             :   ierr = PetscMalloc4(m,&eig,m*m,&Q,m*n,&W,lwork,&work);CHKERRQ(ierr);
     196             : #endif
     197             : 
     198             :   /* compute eigendecomposition */
     199             :   PetscStackCallBLAS("LAPACKlacpy",LAPACKlacpy_("L",&n,&n,A,&ld,Q,&ld));
     200             : #if defined(PETSC_USE_COMPLEX)
     201             :   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
     202             : #else
     203             :   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
     204             : #endif
     205             :   SlepcCheckLapackInfo("syev",info);
     206             : 
     207             :   /* W = f(Lambda)*Q' */
     208             :   for (i=0;i<n;i++) {
     209             :     for (j=0;j<n;j++) W[i+j*ld] = Q[j+i*ld]*PetscAbsScalar(eig[i]);
     210             :   }
     211             :   /* A = Q*W */
     212             :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,W,&ld,&zero,A,&ld));
     213             : #if defined(PETSC_USE_COMPLEX)
     214             :   ierr = PetscFree5(eig,Q,W,work,rwork);CHKERRQ(ierr);
     215             : #else
     216             :   ierr = PetscFree4(eig,Q,W,work);CHKERRQ(ierr);
     217             : #endif
     218             :   PetscFunctionReturn(0);
     219             : #endif
     220             : }
     221             : #endif
     222             : 
     223             : /*
     224             :    Compute the lower Cholesky factor of A
     225             :  */
     226          16 : static PetscErrorCode CholeskyFactor(PetscScalar *A,PetscInt m,PetscInt lda)
     227             : {
     228             : #if defined(PETSC_MISSING_LAPACK_POTRF)
     229             :   PetscFunctionBegin;
     230             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable");
     231             : #else
     232             :   PetscErrorCode ierr;
     233             :   PetscInt       i;
     234             :   PetscScalar    *S;
     235             :   PetscBLASInt   info,n,ld;
     236             : 
     237          16 :   PetscFunctionBegin;
     238          16 :   ierr = PetscBLASIntCast(m,&n);CHKERRQ(ierr);
     239          16 :   ierr = PetscBLASIntCast(lda,&ld);CHKERRQ(ierr);
     240          16 :   ierr = PetscMalloc1(m*m,&S);CHKERRQ(ierr);
     241             : 
     242             :   /* save a copy of matrix in S */
     243         498 :   for (i=0;i<m;i++) {
     244         498 :     ierr = PetscArraycpy(S+i*m,A+i*lda,m);CHKERRQ(ierr);
     245             :   }
     246             : 
     247             :   /* compute upper Cholesky factor in R */
     248          16 :   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,A,&ld,&info));
     249          16 :   ierr = PetscLogFlops((1.0*n*n*n)/3.0);CHKERRQ(ierr);
     250             : 
     251          16 :   if (info) {  /* LAPACKpotrf failed, retry on diagonally perturbed matrix */
     252         498 :     for (i=0;i<m;i++) {
     253         498 :       ierr = PetscArraycpy(A+i*lda,S+i*m,m);CHKERRQ(ierr);
     254         498 :       A[i+i*lda] += 50.0*PETSC_MACHINE_EPSILON;
     255             :     }
     256          16 :     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,A,&ld,&info));
     257          16 :     SlepcCheckLapackInfo("potrf",info);
     258          16 :     ierr = PetscLogFlops((1.0*n*n*n)/3.0);CHKERRQ(ierr);
     259             :   }
     260             : 
     261             :   /* Zero out entries above the diagonal */
     262         482 :   for (i=1;i<m;i++) {
     263         482 :     ierr = PetscArrayzero(A+i*lda,i);CHKERRQ(ierr);
     264             :   }
     265          16 :   ierr = PetscFree(S);CHKERRQ(ierr);
     266          16 :   PetscFunctionReturn(0);
     267             : #endif
     268             : }
     269             : 
     270             : /*
     271             :    LyapunovFact_LAPACK - alternative implementation when SLICOT is not available
     272             : */
     273          16 : static PetscErrorCode LyapunovChol_LAPACK(PetscScalar *H,PetscInt m,PetscInt ldh,PetscScalar *r,PetscScalar *L,PetscInt ldl,PetscReal *res)
     274             : {
     275             : #if defined(PETSC_MISSING_LAPACK_HSEQR) || defined(SLEPC_MISSING_LAPACK_TRSYL)
     276             :   PetscFunctionBegin;
     277             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"HSEQR/TRSYL - Lapack routines are unavailable");
     278             : #else
     279             :   PetscErrorCode ierr;
     280          16 :   PetscBLASInt   ilo=1,lwork,info,n,ld,ld1,ione=1;
     281             :   PetscInt       i,j;
     282             :   PetscReal      scal,beta;
     283          16 :   PetscScalar    *Q,*C,*W,*z,*wr,*work,zero=0.0,done=1.0;
     284             : #if !defined(PETSC_USE_COMPLEX)
     285             :   PetscScalar    *wi;
     286             : #endif
     287             : 
     288          16 :   PetscFunctionBegin;
     289          16 :   ierr = PetscBLASIntCast(ldh,&ld);CHKERRQ(ierr);
     290          16 :   ierr = PetscBLASIntCast(ldl,&ld1);CHKERRQ(ierr);
     291          16 :   ierr = PetscBLASIntCast(m,&n);CHKERRQ(ierr);
     292          16 :   lwork = n;
     293          16 :   C = L;
     294             : 
     295             :   /* compute the (real) Schur form of H */
     296             : #if !defined(PETSC_USE_COMPLEX)
     297          16 :   ierr = PetscMalloc6(m*m,&Q,m*m,&W,m,&z,m,&wr,m,&wi,m,&work);CHKERRQ(ierr);
     298          16 :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,H,&ld,wr,wi,Q,&n,work,&lwork,&info));
     299             : #else
     300             :   ierr = PetscMalloc5(m*m,&Q,m*m,&W,m,&z,m,&wr,m,&work);CHKERRQ(ierr);
     301             :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,H,&ld,wr,Q,&n,work,&lwork,&info));
     302             : #endif
     303          16 :   SlepcCheckLapackInfo("hseqr",info);
     304             : #if defined(PETSC_USE_DEBUG)
     305         498 :   for (i=0;i<m;i++) if (PetscRealPart(wr[i])>0.0) SETERRQ(PETSC_COMM_SELF,1,"Positive eigenvalue found, the coefficient matrix is not stable");
     306             : #endif
     307             : 
     308             :   /* C = z*z', z = Q'*r */
     309          16 :   PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&done,Q,&n,r,&ione,&zero,z,&ione));
     310         498 :   for (i=0;i<m;i++) {
     311       16830 :     for (j=0;j<m;j++) C[i+j*ldl] = -z[i]*PetscConj(z[j]);
     312             :   }
     313             : 
     314             :   /* solve triangular Sylvester equation */
     315          16 :   PetscStackCallBLAS("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,H,&ld,H,&ld,C,&ld1,&scal,&info));
     316          16 :   SlepcCheckLapackInfo("trsyl",info);
     317          16 :   if (scal!=1.0) SETERRQ1(PETSC_COMM_SELF,1,"Current implementation cannot handle scale factor %g",scal);
     318             :   
     319             :   /* back-transform C = Q*C*Q' */
     320          16 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,C,&n,&zero,W,&n));
     321          16 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,C,&ld1));
     322             : 
     323             :   /* resnorm = norm(H(m+1,:)*Y) */
     324          16 :   beta = PetscAbsScalar(H[m+(m-1)*ldh]);
     325          16 :   *res = beta*BLASnrm2_(&n,C+m-1,&n);
     326             : 
     327             : #if 0
     328             :   /* avoid problems due to (small) negative eigenvalues */
     329             :   ierr = AbsEig(C,m);CHKERRQ(ierr);
     330             : #endif
     331             : 
     332             :   /* L = chol(C,'lower') */
     333          16 :   ierr = CholeskyFactor(C,m,ldl);CHKERRQ(ierr);
     334             : 
     335             : #if !defined(PETSC_USE_COMPLEX)
     336          16 :   ierr = PetscFree6(Q,W,z,wr,wi,work);CHKERRQ(ierr);
     337             : #else
     338             :   ierr = PetscFree5(Q,W,z,wr,work);CHKERRQ(ierr);
     339             : #endif
     340          16 :   PetscFunctionReturn(0);
     341             : #endif
     342             : }
     343             : 
     344             : #endif /* SLEPC_HAVE_SLICOT */
     345             : 
     346             : /*@C
     347             :    LMEDenseLyapunovFact - Computes the Cholesky factor of the solution of a
     348             :    dense Lyapunov equation with rank-1 right-hand side.
     349             : 
     350             :    Logically Collective on lme
     351             : 
     352             :    Input Parameters:
     353             : +  lme - linear matrix equation solver context
     354             : .  H   - coefficient matrix
     355             : .  m   - problem size
     356             : .  ldh - leading dimension of H
     357             : .  r   - right-hand side vector
     358             : -  ldl - leading dimension of L
     359             : 
     360             :    Output Parameter:
     361             : +  L   - Cholesky factor of the solution
     362             : -  res - residual norm
     363             : 
     364             :    Note:
     365             :    The Lyapunov equation has the form H*X + X*H' = -r*r', where H represents
     366             :    the leading mxm submatrix of argument H, and the solution X = L*L'.
     367             : 
     368             :    H is assumed to be in upper Hessenberg form, with dimensions (m+1)xm.
     369             :    The last row is used to compute the residual norm, assuming H and r come
     370             :    from the projection onto an Arnoldi basis.
     371             : 
     372             :    Level: developer
     373             : 
     374             : .seealso: LMESolve()
     375             : @*/
     376          16 : PetscErrorCode LMEDenseLyapunovChol(LME lme,PetscScalar *H,PetscInt m,PetscInt ldh,PetscScalar *r,PetscScalar *L,PetscInt ldl,PetscReal *res)
     377             : {
     378             :   PetscErrorCode ierr;
     379             : #if defined(PETSC_USE_INFO)
     380             :   PetscInt       i;
     381             :   PetscScalar    *Hcopy;
     382             :   PetscReal      error;
     383             : #endif
     384             : 
     385          16 :   PetscFunctionBegin;
     386             : #if defined(PETSC_USE_INFO)
     387          16 :   if (PetscLogPrintInfo) {
     388           4 :     ierr = PetscMalloc1(m*m,&Hcopy);CHKERRQ(ierr);
     389         126 :     for (i=0;i<m;i++) {
     390         126 :       ierr = PetscArraycpy(Hcopy+i*m,H+i*ldh,m);CHKERRQ(ierr);
     391             :     }
     392             :   }
     393             : #endif
     394             : 
     395             : #if defined(SLEPC_HAVE_SLICOT)
     396             :   ierr = LyapunovChol_SLICOT(H,m,ldh,r,L,ldl,res);CHKERRQ(ierr);
     397             : #else
     398          16 :   ierr = LyapunovChol_LAPACK(H,m,ldh,r,L,ldl,res);CHKERRQ(ierr);
     399             : #endif
     400             : 
     401             : #if defined(PETSC_USE_INFO)
     402          16 :   if (PetscLogPrintInfo) {
     403           4 :     ierr = LyapunovResidual(Hcopy,m,m,r,L,ldl,&error);CHKERRQ(ierr);
     404           4 :     ierr = PetscInfo1(lme,"Residual norm of dense Lyapunov equation = %g\n",error);CHKERRQ(ierr);
     405           4 :     ierr = PetscFree(Hcopy);CHKERRQ(ierr);
     406             :   }
     407             : #endif
     408          16 :   PetscFunctionReturn(0);
     409             : }
     410             : 

Generated by: LCOV version 1.13