LCOV - code coverage report
Current view: top level - lme/interface - lmesolve.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 122 123 99.2 %
Date: 2019-11-15 09:58:42 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             :    LME routines related to the solution process
      12             : */
      13             : 
      14             : #include <slepc/private/lmeimpl.h>   /*I "slepclme.h" I*/
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : /*@
      18             :    LMESolve - Solves the linear matrix equation.
      19             : 
      20             :    Collective on lme
      21             : 
      22             :    Input Parameters:
      23             : .  lme - linear matrix equation solver context obtained from LMECreate()
      24             : 
      25             :    Options Database Keys:
      26             : +  -lme_view - print information about the solver used
      27             : .  -lme_view_mat binary - save the matrix to the default binary viewer
      28             : .  -lme_view_rhs binary - save right hand side to the default binary viewer
      29             : .  -lme_view_solution binary - save computed solution to the default binary viewer
      30             : -  -lme_converged_reason - print reason for convergence, and number of iterations
      31             : 
      32             :    Notes:
      33             :    The matrix coefficients are specified with LMESetCoefficients().
      34             :    The right-hand side is specified with LMESetRHS().
      35             :    The placeholder for the solution is specified with LMESetSolution().
      36             : 
      37             :    Level: beginner
      38             : 
      39             : .seealso: LMECreate(), LMESetUp(), LMEDestroy(), LMESetTolerances(), LMESetCoefficients(), LMESetRHS(), LMESetSolution()
      40             : @*/
      41           5 : PetscErrorCode LMESolve(LME lme)
      42             : {
      43             :   PetscErrorCode ierr;
      44             : 
      45           5 :   PetscFunctionBegin;
      46           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
      47             : 
      48             :   /* call setup */
      49           5 :   ierr = LMESetUp(lme);CHKERRQ(ierr);
      50           5 :   lme->its    = 0;
      51           5 :   lme->errest = 0.0;
      52             : 
      53           5 :   ierr = LMEViewFromOptions(lme,NULL,"-lme_view_pre");CHKERRQ(ierr);
      54             : 
      55             :   /* call solver */
      56           5 :   if (!lme->ops->solve[lme->problem_type]) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"The specified solver does not support equation type %s",LMEProblemTypes[lme->problem_type]);
      57           5 :   ierr = PetscLogEventBegin(LME_Solve,lme,0,0,0);CHKERRQ(ierr);
      58           5 :   ierr = (*lme->ops->solve[lme->problem_type])(lme);CHKERRQ(ierr);
      59           5 :   ierr = PetscLogEventEnd(LME_Solve,lme,0,0,0);CHKERRQ(ierr);
      60             : 
      61           5 :   if (!lme->reason) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
      62             : 
      63           5 :   if (lme->errorifnotconverged && lme->reason < 0) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_NOT_CONVERGED,"LMESolve has not converged");
      64             : 
      65             :   /* various viewers */
      66           5 :   ierr = LMEViewFromOptions(lme,NULL,"-lme_view");CHKERRQ(ierr);
      67           5 :   ierr = LMEReasonViewFromOptions(lme);CHKERRQ(ierr);
      68           5 :   ierr = MatViewFromOptions(lme->A,(PetscObject)lme,"-lme_view_mat");CHKERRQ(ierr);
      69           5 :   ierr = MatViewFromOptions(lme->C,(PetscObject)lme,"-lme_view_rhs");CHKERRQ(ierr);
      70           5 :   ierr = MatViewFromOptions(lme->X,(PetscObject)lme,"-lme_view_solution");CHKERRQ(ierr);
      71           5 :   PetscFunctionReturn(0);
      72             : }
      73             : 
      74             : /*@
      75             :    LMEGetIterationNumber - Gets the current iteration number. If the
      76             :    call to LMESolve() is complete, then it returns the number of iterations
      77             :    carried out by the solution method.
      78             : 
      79             :    Not Collective
      80             : 
      81             :    Input Parameter:
      82             : .  lme - the linear matrix equation solver context
      83             : 
      84             :    Output Parameter:
      85             : .  its - number of iterations
      86             : 
      87             :    Note:
      88             :    During the i-th iteration this call returns i-1. If LMESolve() is
      89             :    complete, then parameter "its" contains either the iteration number at
      90             :    which convergence was successfully reached, or failure was detected.
      91             :    Call LMEGetConvergedReason() to determine if the solver converged or
      92             :    failed and why.
      93             : 
      94             :    Level: intermediate
      95             : 
      96             : .seealso: LMEGetConvergedReason(), LMESetTolerances()
      97             : @*/
      98           2 : PetscErrorCode LMEGetIterationNumber(LME lme,PetscInt *its)
      99             : {
     100           2 :   PetscFunctionBegin;
     101           2 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     102           2 :   PetscValidIntPointer(its,2);
     103           2 :   *its = lme->its;
     104           2 :   PetscFunctionReturn(0);
     105             : }
     106             : 
     107             : /*@
     108             :    LMEGetConvergedReason - Gets the reason why the LMESolve() iteration was
     109             :    stopped.
     110             : 
     111             :    Not Collective
     112             : 
     113             :    Input Parameter:
     114             : .  lme - the linear matrix equation solver context
     115             : 
     116             :    Output Parameter:
     117             : .  reason - negative value indicates diverged, positive value converged
     118             : 
     119             :    Notes:
     120             : 
     121             :    Possible values for reason are
     122             : +  LME_CONVERGED_TOL - converged up to tolerance
     123             : .  LME_DIVERGED_ITS - required more than max_it iterations to reach convergence
     124             : -  LME_DIVERGED_BREAKDOWN - generic breakdown in method
     125             : 
     126             :    Can only be called after the call to LMESolve() is complete.
     127             : 
     128             :    Level: intermediate
     129             : 
     130             : .seealso: LMESetTolerances(), LMESolve(), LMEConvergedReason, LMESetErrorIfNotConverged()
     131             : @*/
     132           2 : PetscErrorCode LMEGetConvergedReason(LME lme,LMEConvergedReason *reason)
     133             : {
     134           2 :   PetscFunctionBegin;
     135           2 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     136           2 :   PetscValidIntPointer(reason,2);
     137           2 :   *reason = lme->reason;
     138           2 :   PetscFunctionReturn(0);
     139             : }
     140             : 
     141             : /*@
     142             :    LMEGetErrorEstimate - Returns the error estimate obtained during solve.
     143             : 
     144             :    Not Collective
     145             : 
     146             :    Input Parameter:
     147             : .  lme - linear matrix equation solver context
     148             : 
     149             :    Output Parameter:
     150             : .  errest - the error estimate
     151             : 
     152             :    Notes:
     153             :    This is the error estimated internally by the solver. The actual
     154             :    error bound can be computed with LMEComputeError(). Note that some
     155             :    solvers may not be able to provide an error estimate.
     156             : 
     157             :    Level: advanced
     158             : 
     159             : .seealso: LMEComputeError()
     160             : @*/
     161           5 : PetscErrorCode LMEGetErrorEstimate(LME lme,PetscReal *errest)
     162             : {
     163           5 :   PetscFunctionBegin;
     164           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     165           5 :   PetscValidRealPointer(errest,2);
     166           5 :   *errest = lme->errest;
     167           5 :   PetscFunctionReturn(0);
     168             : }
     169             : 
     170             : /*
     171             :    LMEComputeResidualNorm_Lyapunov - Computes the Frobenius norm of the residual matrix
     172             :    associated with the Lyapunov equation.
     173             : */
     174           5 : PetscErrorCode LMEComputeResidualNorm_Lyapunov(LME lme,PetscReal *norm)
     175             : {
     176             :   PetscErrorCode    ierr;
     177             :   PetscInt          j,n,N,k,l;
     178             :   PetscBLASInt      n_,N_,k_,l_;
     179           5 :   PetscScalar       *Rarray,alpha=1.0,beta=0.0;
     180             :   const PetscScalar *A,*B;
     181             :   BV                W,AX,X1,C1;
     182             :   Mat               R,X1m,C1m;
     183             :   Vec               v,w;
     184             :   VecScatter        vscat;
     185             : 
     186           5 :   PetscFunctionBegin;
     187           5 :   ierr = MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL);CHKERRQ(ierr);
     188           5 :   ierr = MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL);CHKERRQ(ierr);
     189           5 :   ierr = BVCreateFromMat(C1m,&C1);CHKERRQ(ierr);
     190           5 :   ierr = BVSetFromOptions(C1);CHKERRQ(ierr);
     191           5 :   ierr = BVCreateFromMat(X1m,&X1);CHKERRQ(ierr);
     192           5 :   ierr = BVSetFromOptions(X1);CHKERRQ(ierr);
     193           5 :   ierr = BVGetSizes(X1,&n,&N,&k);CHKERRQ(ierr);
     194           5 :   ierr = BVGetSizes(C1,NULL,NULL,&l);CHKERRQ(ierr);
     195           5 :   ierr = PetscBLASIntCast(n,&n_);CHKERRQ(ierr);
     196           5 :   ierr = PetscBLASIntCast(N,&N_);CHKERRQ(ierr);
     197           5 :   ierr = PetscBLASIntCast(k,&k_);CHKERRQ(ierr);
     198           5 :   ierr = PetscBLASIntCast(l,&l_);CHKERRQ(ierr);
     199             : 
     200             :   /* create W to store a redundant copy of a BV in each process */
     201           5 :   ierr = BVCreate(PETSC_COMM_SELF,&W);CHKERRQ(ierr);
     202           5 :   ierr = BVSetSizes(W,N,N,k);CHKERRQ(ierr);
     203           5 :   ierr = BVSetFromOptions(W);CHKERRQ(ierr);
     204           5 :   ierr = BVGetColumn(X1,0,&v);CHKERRQ(ierr);
     205           5 :   ierr = VecScatterCreateToAll(v,&vscat,NULL);CHKERRQ(ierr);
     206           5 :   ierr = BVRestoreColumn(X1,0,&v);CHKERRQ(ierr);
     207             : 
     208             :   /* create AX to hold the product A*X1 */
     209           5 :   ierr = BVDuplicate(X1,&AX);CHKERRQ(ierr);
     210           5 :   ierr = BVMatMult(X1,lme->A,AX);CHKERRQ(ierr);
     211             : 
     212             :   /* create dense matrix to hold the residual R=C1*C1'+AX*X1'+X1*AX' */
     213           5 :   ierr = MatCreateDense(PetscObjectComm((PetscObject)lme),n,n,N,N,NULL,&R);CHKERRQ(ierr);
     214             : 
     215             :   /* R=C1*C1' */
     216           5 :   ierr = MatDenseGetArray(R,&Rarray);CHKERRQ(ierr);
     217          10 :   for (j=0;j<l;j++) {
     218          10 :     ierr = BVGetColumn(C1,j,&v);CHKERRQ(ierr);
     219          10 :     ierr = BVGetColumn(W,j,&w);CHKERRQ(ierr);
     220          10 :     ierr = VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     221          10 :     ierr = VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     222          10 :     ierr = BVRestoreColumn(C1,j,&v);CHKERRQ(ierr);
     223          10 :     ierr = BVRestoreColumn(W,j,&w);CHKERRQ(ierr);
     224             :   }
     225           5 :   if (n) {
     226           5 :     ierr = BVGetArrayRead(C1,&A);CHKERRQ(ierr);
     227           5 :     ierr = BVGetArrayRead(W,&B);CHKERRQ(ierr);
     228           5 :     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&l_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
     229           5 :     ierr = BVRestoreArrayRead(C1,&A);CHKERRQ(ierr);
     230           5 :     ierr = BVRestoreArrayRead(W,&B);CHKERRQ(ierr);
     231             :   }
     232           5 :   beta = 1.0;
     233             : 
     234             :   /* R+=AX*X1' */
     235         357 :   for (j=0;j<k;j++) {
     236         352 :     ierr = BVGetColumn(X1,j,&v);CHKERRQ(ierr);
     237         352 :     ierr = BVGetColumn(W,j,&w);CHKERRQ(ierr);
     238         352 :     ierr = VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     239         352 :     ierr = VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     240         352 :     ierr = BVRestoreColumn(X1,j,&v);CHKERRQ(ierr);
     241         352 :     ierr = BVRestoreColumn(W,j,&w);CHKERRQ(ierr);
     242             :   }
     243           5 :   if (n) {
     244           5 :     ierr = BVGetArrayRead(AX,&A);CHKERRQ(ierr);
     245           5 :     ierr = BVGetArrayRead(W,&B);CHKERRQ(ierr);
     246           5 :     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
     247           5 :     ierr = BVRestoreArrayRead(AX,&A);CHKERRQ(ierr);
     248           5 :     ierr = BVRestoreArrayRead(W,&B);CHKERRQ(ierr);
     249             :   }
     250             : 
     251             :   /* R+=X1*AX' */
     252         352 :   for (j=0;j<k;j++) {
     253         352 :     ierr = BVGetColumn(AX,j,&v);CHKERRQ(ierr);
     254         352 :     ierr = BVGetColumn(W,j,&w);CHKERRQ(ierr);
     255         352 :     ierr = VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     256         352 :     ierr = VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     257         352 :     ierr = BVRestoreColumn(AX,j,&v);CHKERRQ(ierr);
     258         352 :     ierr = BVRestoreColumn(W,j,&w);CHKERRQ(ierr);
     259             :   }
     260           5 :   if (n) {
     261           5 :     ierr = BVGetArrayRead(X1,&A);CHKERRQ(ierr);
     262           5 :     ierr = BVGetArrayRead(W,&B);CHKERRQ(ierr);
     263           5 :     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
     264           5 :     ierr = BVRestoreArrayRead(X1,&A);CHKERRQ(ierr);
     265           5 :     ierr = BVRestoreArrayRead(W,&B);CHKERRQ(ierr);
     266             :   }
     267           5 :   ierr = MatDenseRestoreArray(R,&Rarray);CHKERRQ(ierr);
     268             : 
     269             :   /* compute ||R||_F */
     270           5 :   ierr = MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     271           5 :   ierr = MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     272           5 :   ierr = MatNorm(R,NORM_FROBENIUS,norm);CHKERRQ(ierr);
     273             : 
     274           5 :   ierr = BVDestroy(&W);CHKERRQ(ierr);
     275           5 :   ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr);
     276           5 :   ierr = BVDestroy(&AX);CHKERRQ(ierr);
     277           5 :   ierr = MatDestroy(&R);CHKERRQ(ierr);
     278           5 :   ierr = BVDestroy(&C1);CHKERRQ(ierr);
     279           5 :   ierr = BVDestroy(&X1);CHKERRQ(ierr);
     280           5 :   PetscFunctionReturn(0);
     281             : }
     282             : 
     283             : /*@
     284             :    LMEComputeError - Computes the error (based on the residual norm) associated
     285             :    with the last equation solved.
     286             : 
     287             :    Collective on lme
     288             : 
     289             :    Input Parameter:
     290             : .  lme  - the linear matrix equation solver context
     291             : 
     292             :    Output Parameter:
     293             : .  error - the error
     294             : 
     295             :    Notes:
     296             :    This function is not scalable (in terms of memory or parallel communication),
     297             :    so it should not be called except in the case of small problem size. For
     298             :    large equations, use LMEGetErrorEstimate().
     299             : 
     300             :    Level: advanced
     301             : 
     302             : .seealso: LMESolve(), LMEGetErrorEstimate()
     303             : @*/
     304           5 : PetscErrorCode LMEComputeError(LME lme,PetscReal *error)
     305             : {
     306             :   PetscErrorCode ierr;
     307             : 
     308           5 :   PetscFunctionBegin;
     309           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     310           5 :   PetscValidRealPointer(error,2);
     311             : 
     312           5 :   ierr = PetscLogEventBegin(LME_ComputeError,lme,0,0,0);CHKERRQ(ierr);
     313             :   /* compute residual norm */
     314           5 :   switch (lme->problem_type) {
     315             :     case LME_LYAPUNOV:
     316           5 :       ierr = LMEComputeResidualNorm_Lyapunov(lme,error);CHKERRQ(ierr);
     317             :       break;
     318             :     default:
     319           0 :       SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"Not implemented for equation type %s",LMEProblemTypes[lme->problem_type]);
     320             :   }
     321             : 
     322             :   /* compute error */
     323             :   /* currently we only support absolute error, so just return the norm */
     324           5 :   ierr = PetscLogEventEnd(LME_ComputeError,lme,0,0,0);CHKERRQ(ierr);
     325           5 :   PetscFunctionReturn(0);
     326             : }
     327             : 

Generated by: LCOV version 1.13