LCOV - code coverage report
Current view: top level - lme/interface - lmesetup.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 112 127 88.2 %
Date: 2019-11-15 09:58:42 Functions: 10 10 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 problem setup
      12             : */
      13             : 
      14             : #include <slepc/private/lmeimpl.h>       /*I "slepclme.h" I*/
      15             : 
      16           5 : PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Lyapunov(LME lme)
      17             : {
      18             :   PetscErrorCode ierr;
      19             :   Mat            C1,C2,X1,X2;
      20             :   Vec            dc,dx;
      21             : 
      22           5 :   PetscFunctionBegin;
      23           5 :   ierr = MatLRCGetMats(lme->C,NULL,&C1,&dc,&C2);CHKERRQ(ierr);
      24           5 :   if (C1!=C2) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric right-hand side C");
      25           5 :   if (dc) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently require positive-definite right-hand side C");
      26           5 :   if (lme->X) {
      27           1 :     ierr = MatLRCGetMats(lme->X,NULL,&X1,&dx,&X2);CHKERRQ(ierr);
      28           1 :     if (X1!=X2) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric solution X");
      29           1 :     if (dx) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently assume a positive-definite solution X");
      30             :   }
      31           5 :   PetscFunctionReturn(0);
      32             : }
      33             : 
      34             : /*@
      35             :    LMESetUp - Sets up all the internal data structures necessary for the
      36             :    execution of the linear matrix equation solver.
      37             : 
      38             :    Collective on lme
      39             : 
      40             :    Input Parameter:
      41             : .  lme   - linear matrix equation solver context
      42             : 
      43             :    Notes:
      44             :    This function need not be called explicitly in most cases, since LMESolve()
      45             :    calls it. It can be useful when one wants to measure the set-up time
      46             :    separately from the solve time.
      47             : 
      48             :    Level: developer
      49             : 
      50             : .seealso: LMECreate(), LMESolve(), LMEDestroy()
      51             : @*/
      52           5 : PetscErrorCode LMESetUp(LME lme)
      53             : {
      54             :   PetscErrorCode ierr;
      55             :   PetscInt       N;
      56             : 
      57           5 :   PetscFunctionBegin;
      58           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
      59             : 
      60             :   /* reset the convergence flag from the previous solves */
      61           5 :   lme->reason = LME_CONVERGED_ITERATING;
      62             : 
      63           5 :   if (lme->setupcalled) PetscFunctionReturn(0);
      64           5 :   ierr = PetscLogEventBegin(LME_SetUp,lme,0,0,0);CHKERRQ(ierr);
      65             : 
      66             :   /* Set default solver type (LMESetFromOptions was not called) */
      67           5 :   if (!((PetscObject)lme)->type_name) {
      68           0 :     ierr = LMESetType(lme,LMEKRYLOV);CHKERRQ(ierr);
      69             :   }
      70             : 
      71             :   /* Check problem dimensions */
      72           5 :   if (!lme->A) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
      73           5 :   ierr = MatGetSize(lme->A,&N,NULL);CHKERRQ(ierr);
      74           5 :   if (lme->ncv > N) lme->ncv = N;
      75             : 
      76             :   /* setup options for the particular equation type */
      77           5 :   switch (lme->problem_type) {
      78             :     case LME_LYAPUNOV:
      79           5 :       ierr = LMESetUp_Lyapunov(lme);CHKERRQ(ierr);
      80             :       break;
      81             :     case LME_SYLVESTER:
      82           0 :       LMECheckCoeff(lme,lme->B,"B","Sylvester");
      83             :       break;
      84             :     case LME_GEN_LYAPUNOV:
      85           0 :       LMECheckCoeff(lme,lme->D,"D","Generalized Lyapunov");
      86             :       break;
      87             :     case LME_GEN_SYLVESTER:
      88           0 :       LMECheckCoeff(lme,lme->B,"B","Generalized Sylvester");
      89           0 :       LMECheckCoeff(lme,lme->D,"D","Generalized Sylvester");
      90           0 :       LMECheckCoeff(lme,lme->E,"E","Generalized Sylvester");
      91             :       break;
      92             :     case LME_DT_LYAPUNOV:
      93             :       break;
      94             :     case LME_STEIN:
      95           0 :       LMECheckCoeff(lme,lme->D,"D","Stein");
      96             :       break;
      97             :   }
      98           5 :   if (lme->problem_type!=LME_LYAPUNOV) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
      99             : 
     100             :   /* call specific solver setup */
     101           5 :   ierr = (*lme->ops->setup)(lme);CHKERRQ(ierr);
     102             : 
     103             :   /* set tolerance if not yet set */
     104           5 :   if (lme->tol==PETSC_DEFAULT) lme->tol = SLEPC_DEFAULT_TOL;
     105             : 
     106           5 :   ierr = PetscLogEventEnd(LME_SetUp,lme,0,0,0);CHKERRQ(ierr);
     107           5 :   lme->setupcalled = 1;
     108           5 :   PetscFunctionReturn(0);
     109             : }
     110             : 
     111           5 : PETSC_STATIC_INLINE PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
     112             : {
     113             :   PetscErrorCode ierr;
     114             :   PetscInt       m,n;
     115             : 
     116           5 :   PetscFunctionBegin;
     117           5 :   ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
     118           5 :   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
     119           5 :   if (!lme->setupcalled) { ierr = MatDestroy(lmeA);CHKERRQ(ierr); }
     120           5 :   ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
     121           5 :   *lmeA = A;
     122           5 :   PetscFunctionReturn(0);
     123             : }
     124             : 
     125             : /*@
     126             :    LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
     127             :    equation to be solved.
     128             : 
     129             :    Collective on lme
     130             : 
     131             :    Input Parameters:
     132             : +  lme - the matrix function context
     133             : .  A   - first coefficient matrix
     134             : .  B   - second coefficient matrix
     135             : .  D   - third coefficient matrix
     136             : -  E   - fourth coefficient matrix
     137             : 
     138             :    Notes:
     139             :    The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is not
     140             :    provided here but with LMESetRHS(). Not all four matrices must be passed, some
     141             :    can be NULL instead, see LMESetProblemType() for details.
     142             : 
     143             :    It must be called before LMESetUp(). If it is called again after LMESetUp() then
     144             :    the LME object is reset.
     145             : 
     146             :    In order to delete a previously set matrix, pass a NULL in the corresponding
     147             :    argument.
     148             : 
     149             :    Level: beginner
     150             : 
     151             : .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
     152             : @*/
     153           5 : PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
     154             : {
     155             :   PetscErrorCode ierr;
     156             : 
     157           5 :   PetscFunctionBegin;
     158           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     159           5 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     160           5 :   PetscCheckSameComm(lme,1,A,2);
     161           5 :   if (B) {
     162           0 :     PetscValidHeaderSpecific(B,MAT_CLASSID,3);
     163           0 :     PetscCheckSameComm(lme,1,B,3);
     164             :   }
     165           5 :   if (D) {
     166           0 :     PetscValidHeaderSpecific(D,MAT_CLASSID,4);
     167           0 :     PetscCheckSameComm(lme,1,D,4);
     168             :   }
     169           5 :   if (E) {
     170           0 :     PetscValidHeaderSpecific(E,MAT_CLASSID,5);
     171           0 :     PetscCheckSameComm(lme,1,E,5);
     172             :   }
     173             : 
     174           5 :   if (lme->setupcalled) { ierr = LMEReset(lme);CHKERRQ(ierr); }
     175             : 
     176           5 :   ierr = LMESetCoefficients_Private(lme,A,&lme->A);CHKERRQ(ierr);
     177           5 :   if (B) { ierr = LMESetCoefficients_Private(lme,B,&lme->B);CHKERRQ(ierr); }
     178           5 :   else if (!lme->setupcalled) { ierr = MatDestroy(&lme->B);CHKERRQ(ierr); }
     179           5 :   if (D) { ierr = LMESetCoefficients_Private(lme,D,&lme->D);CHKERRQ(ierr); }
     180           5 :   else if (!lme->setupcalled) { ierr = MatDestroy(&lme->D);CHKERRQ(ierr); }
     181           5 :   if (E) { ierr = LMESetCoefficients_Private(lme,E,&lme->E);CHKERRQ(ierr); }
     182           5 :   else if (!lme->setupcalled) { ierr = MatDestroy(&lme->E);CHKERRQ(ierr); }
     183             : 
     184           5 :   lme->setupcalled = 0;
     185           5 :   PetscFunctionReturn(0);
     186             : }
     187             : 
     188             : /*@
     189             :    LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.
     190             : 
     191             :    Collective on lme
     192             : 
     193             :    Input Parameter:
     194             : .  lme - the LME context
     195             : 
     196             :    Output Parameters:
     197             : +  A   - first coefficient matrix
     198             : .  B   - second coefficient matrix
     199             : .  D   - third coefficient matrix
     200             : -  E   - fourth coefficient matrix
     201             : 
     202             :    Level: intermediate
     203             : 
     204             : .seealso: LMESolve(), LMESetCoefficients()
     205             : @*/
     206           3 : PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
     207             : {
     208           3 :   PetscFunctionBegin;
     209           3 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     210           3 :   if (A) *A = lme->A;
     211           3 :   if (B) *B = lme->B;
     212           3 :   if (D) *D = lme->D;
     213           3 :   if (E) *E = lme->E;
     214           3 :   PetscFunctionReturn(0);
     215             : }
     216             : 
     217             : /*@
     218             :    LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
     219             :    matrix.
     220             : 
     221             :    Collective on lme
     222             : 
     223             :    Input Parameters:
     224             : +  lme - the matrix function context
     225             : -  C   - the right-hand side matrix
     226             : 
     227             :    Notes:
     228             :    The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is
     229             :    given with this function. C must be a low-rank matrix of type MATLRC, that is,
     230             :    C = U*D*V' where D is diagonal of order k, and U, V are dense tall-skinny
     231             :    matrices with k columns. No sparse matrix must be provided when creating the
     232             :    MATLRC matrix.
     233             : 
     234             :    In equation types that require C to be symmetric, such as Lyapunov, C must be
     235             :    created with V=U (or V=NULL).
     236             : 
     237             :    Level: beginner
     238             : 
     239             : .seealso: LMESetSolution(), LMESetProblemType()
     240             : @*/
     241           5 : PetscErrorCode LMESetRHS(LME lme,Mat C)
     242             : {
     243             :   PetscErrorCode ierr;
     244             :   PetscBool      match;
     245             :   Mat            A;
     246             : 
     247           5 :   PetscFunctionBegin;
     248           5 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     249           5 :   PetscValidHeaderSpecific(C,MAT_CLASSID,2);
     250           5 :   PetscCheckSameComm(lme,1,C,2);
     251             : 
     252           5 :   ierr = PetscObjectTypeCompare((PetscObject)C,MATLRC,&match);CHKERRQ(ierr);
     253           5 :   if (!match) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
     254           5 :   ierr = MatLRCGetMats(C,&A,NULL,NULL,NULL);CHKERRQ(ierr);
     255           5 :   if (A) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
     256             : 
     257           5 :   ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
     258           5 :   ierr = MatDestroy(&lme->C);CHKERRQ(ierr);
     259           5 :   lme->C = C;
     260           5 :   PetscFunctionReturn(0);
     261             : }
     262             : 
     263             : /*@
     264             :    LMEGetRHS - Gets the right-hand side of the matrix equation.
     265             : 
     266             :    Collective on lme
     267             : 
     268             :    Input Parameter:
     269             : .  lme - the LME context
     270             : 
     271             :    Output Parameters:
     272             : .  C   - the low-rank matrix
     273             : 
     274             :    Level: intermediate
     275             : 
     276             : .seealso: LMESolve(), LMESetRHS()
     277             : @*/
     278           3 : PetscErrorCode LMEGetRHS(LME lme,Mat *C)
     279             : {
     280           3 :   PetscFunctionBegin;
     281           3 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     282           3 :   PetscValidPointer(C,2);
     283           3 :   *C = lme->C;
     284           3 :   PetscFunctionReturn(0);
     285             : }
     286             : 
     287             : /*@
     288             :    LMESetSolution - Sets the placeholder for the solution of the matrix
     289             :    equation, as a low-rank matrix.
     290             : 
     291             :    Collective on lme
     292             : 
     293             :    Input Parameters:
     294             : +  lme - the matrix function context
     295             : -  X   - the solution matrix
     296             : 
     297             :    Notes:
     298             :    The matrix equation takes the general form A*X*E+D*X*B=C, where the solution
     299             :    matrix is of low rank and is written in factored form X = U*D*V'. This function
     300             :    provides a Mat object of type MATLRC that stores U, V and (optionally) D.
     301             :    These factors will be computed during LMESolve().
     302             : 
     303             :    In equation types whose solution X is symmetric, such as Lyapunov, X must be
     304             :    created with V=U (or V=NULL).
     305             : 
     306             :    If the user provides X with this function, then the solver will
     307             :    return a solution with rank at most the number of columns of U. Alternatively,
     308             :    it is possible to let the solver choose the rank of the solution, by
     309             :    setting X to NULL and then calling LMEGetSolution() after LMESolve().
     310             : 
     311             :    Level: intermediate
     312             : 
     313             : .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
     314             : @*/
     315           1 : PetscErrorCode LMESetSolution(LME lme,Mat X)
     316             : {
     317             :   PetscErrorCode ierr;
     318             :   PetscBool      match;
     319             :   Mat            A;
     320             : 
     321           1 :   PetscFunctionBegin;
     322           1 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     323           1 :   if (X) {
     324           1 :     PetscValidHeaderSpecific(X,MAT_CLASSID,2);
     325           1 :     PetscCheckSameComm(lme,1,X,2);
     326           1 :     ierr = PetscObjectTypeCompare((PetscObject)X,MATLRC,&match);CHKERRQ(ierr);
     327           1 :     if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
     328           1 :     ierr = MatLRCGetMats(X,&A,NULL,NULL,NULL);CHKERRQ(ierr);
     329           1 :     if (A) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
     330           1 :     ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr);
     331             :   }
     332           1 :   ierr = MatDestroy(&lme->X);CHKERRQ(ierr);
     333           1 :   lme->X = X;
     334           1 :   PetscFunctionReturn(0);
     335             : }
     336             : 
     337             : /*@
     338             :    LMEGetSolution - Gets the solution of the matrix equation.
     339             : 
     340             :    Collective on lme
     341             : 
     342             :    Input Parameter:
     343             : .  lme - the LME context
     344             : 
     345             :    Output Parameters:
     346             : .  X   - the low-rank matrix
     347             : 
     348             :    Level: intermediate
     349             : 
     350             : .seealso: LMESolve(), LMESetSolution()
     351             : @*/
     352           1 : PetscErrorCode LMEGetSolution(LME lme,Mat *X)
     353             : {
     354           1 :   PetscFunctionBegin;
     355           1 :   PetscValidHeaderSpecific(lme,LME_CLASSID,1);
     356           1 :   PetscValidPointer(X,2);
     357           1 :   *X = lme->X;
     358           1 :   PetscFunctionReturn(0);
     359             : }
     360             : 
     361             : /*@
     362             :    LMEAllocateSolution - Allocate memory storage for common variables such
     363             :    as the basis vectors.
     364             : 
     365             :    Collective on lme
     366             : 
     367             :    Input Parameters:
     368             : +  lme   - linear matrix equation solver context
     369             : -  extra - number of additional positions, used for methods that require a
     370             :            working basis slightly larger than ncv
     371             : 
     372             :    Developers Note:
     373             :    This is SLEPC_EXTERN because it may be required by user plugin LME
     374             :    implementations.
     375             : 
     376             :    Level: developer
     377             : @*/
     378           5 : PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
     379             : {
     380             :   PetscErrorCode ierr;
     381             :   PetscInt       oldsize,requested;
     382             :   Vec            t;
     383             : 
     384           5 :   PetscFunctionBegin;
     385           5 :   requested = lme->ncv + extra;
     386             : 
     387             :   /* oldsize is zero if this is the first time setup is called */
     388           5 :   ierr = BVGetSizes(lme->V,NULL,NULL,&oldsize);CHKERRQ(ierr);
     389             : 
     390             :   /* allocate basis vectors */
     391           5 :   if (!lme->V) { ierr = LMEGetBV(lme,&lme->V);CHKERRQ(ierr); }
     392           5 :   if (!oldsize) {
     393           5 :     if (!((PetscObject)(lme->V))->type_name) {
     394           0 :       ierr = BVSetType(lme->V,BVSVEC);CHKERRQ(ierr);
     395             :     }
     396           5 :     ierr = MatCreateVecsEmpty(lme->A,&t,NULL);CHKERRQ(ierr);
     397           5 :     ierr = BVSetSizesFromVec(lme->V,t,requested);CHKERRQ(ierr);
     398           5 :     ierr = VecDestroy(&t);CHKERRQ(ierr);
     399             :   } else {
     400           0 :     ierr = BVResize(lme->V,requested,PETSC_FALSE);CHKERRQ(ierr);
     401             :   }
     402           5 :   PetscFunctionReturn(0);
     403             : }
     404             : 

Generated by: LCOV version 1.13