Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, 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             :    BV bi-orthogonalization routines
      12             : */
      13             : 
      14             : #include <slepc/private/bvimpl.h>          /*I   "slepcbv.h"   I*/
      15             : 
      16             : /*
      17             :    BVBiorthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt bi-orthogonalization
      18             : */
      19         112 : static PetscErrorCode BVBiorthogonalizeMGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
      20             : {
      21         112 :   PetscInt       i;
      22         112 :   PetscScalar    dot;
      23         112 :   Vec            vi,wi;
      24             : 
      25         112 :   PetscFunctionBegin;
      26         448 :   for (i=-V->nc;i<V->k;i++) {
      27         336 :     PetscCall(BVGetColumn(W,i,&wi));
      28             :     /* h_i = (v, w_i) */
      29         336 :     PetscCall(VecDot(v,wi,&dot));
      30         336 :     PetscCall(BVRestoreColumn(W,i,&wi));
      31             :     /* v <- v - h_i v_i */
      32         336 :     PetscCall(BV_SetValue(V,i,0,c,dot));
      33         336 :     PetscCall(BVGetColumn(V,i,&vi));
      34         336 :     PetscCall(VecAXPY(v,-dot,vi));
      35         336 :     PetscCall(BVRestoreColumn(V,i,&vi));
      36             :   }
      37         112 :   PetscCall(BV_AddCoefficients(V,V->k,h,c));
      38         112 :   PetscFunctionReturn(PETSC_SUCCESS);
      39             : }
      40             : 
      41             : /*
      42             :    BVBiorthogonalizeCGS1 - Compute one step of CGS bi-orthogonalization: v = (I-V*W')*v
      43             : */
      44        1756 : static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
      45             : {
      46        1756 :   PetscFunctionBegin;
      47             :   /* h = W'*v */
      48        1756 :   PetscCall(BVDotVec(W,v,c));
      49             : 
      50             :   /* v = v - V h */
      51        1756 :   PetscCall(BVMultVec(V,-1.0,1.0,v,c));
      52             : 
      53        1756 :   PetscCall(BV_AddCoefficients(V,V->k,h,c));
      54        1756 :   PetscFunctionReturn(PETSC_SUCCESS);
      55             : }
      56             : 
      57             : #define BVBiorthogonalizeGS1(a,b,c,d,e) ((V->orthog_type==BV_ORTHOG_MGS)?BVBiorthogonalizeMGS1:BVBiorthogonalizeCGS1)(a,b,c,d,e)
      58             : 
      59             : /*
      60             :    BVBiorthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
      61             : 
      62             :    V, W - the two basis vectors objects
      63             :    v    - the vector to bi-orthogonalize
      64             : */
      65         934 : static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
      66             : {
      67         934 :   PetscScalar    *h,*c;
      68             : 
      69         934 :   PetscFunctionBegin;
      70         934 :   h = V->h;
      71         934 :   c = V->c;
      72         934 :   PetscCall(BV_CleanCoefficients(V,V->k,h));
      73        1812 :   PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
      74        1812 :   if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
      75         934 :   PetscFunctionReturn(PETSC_SUCCESS);
      76             : }
      77             : 
      78             : /*@
      79             :    BVBiorthogonalizeColumn - Bi-orthogonalize a column of two BV objects.
      80             : 
      81             :    Collective
      82             : 
      83             :    Input Parameters:
      84             : +  V - first basis vectors context
      85             : .  W - second basis vectors context
      86             : -  j - index of column to be bi-orthonormalized
      87             : 
      88             :    Notes:
      89             :    This function bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
      90             :    and V[0..j-1], respectively, so that W[0..j]'*V[0..j] = diagonal.
      91             : 
      92             :    Level: advanced
      93             : 
      94             : .seealso: BVOrthogonalizeColumn(), BVBiorthonormalizeColumn()
      95             : @*/
      96         394 : PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
      97             : {
      98         394 :   PetscInt       ksavev,lsavev,ksavew,lsavew;
      99         394 :   Vec            y,z;
     100             : 
     101         394 :   PetscFunctionBegin;
     102         394 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     103         394 :   PetscValidHeaderSpecific(W,BV_CLASSID,2);
     104        1182 :   PetscValidLogicalCollectiveInt(V,j,3);
     105         394 :   PetscValidType(V,1);
     106         394 :   BVCheckSizes(V,1);
     107         394 :   PetscValidType(W,2);
     108         394 :   BVCheckSizes(W,2);
     109         394 :   PetscCheckSameTypeAndComm(V,1,W,2);
     110         394 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     111         394 :   PetscCheck(j<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but V only has %" PetscInt_FMT " columns",j,V->m);
     112         394 :   PetscCheck(j<W->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but W only has %" PetscInt_FMT " columns",j,W->m);
     113         394 :   PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
     114         394 :   PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
     115         394 :   PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
     116         394 :   PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
     117             : 
     118             :   /* bi-orthogonalize */
     119         394 :   PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
     120         394 :   ksavev = V->k;
     121         394 :   lsavev = V->l;
     122         394 :   ksavew = W->k;
     123         394 :   lsavew = W->l;
     124         394 :   V->k = j;
     125         394 :   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
     126         394 :   W->k = j;
     127         394 :   W->l = -W->nc;
     128         394 :   PetscCall(BV_AllocateCoeffs(V));
     129         394 :   PetscCall(BV_AllocateCoeffs(W));
     130         394 :   PetscCall(BVGetColumn(V,j,&y));
     131         394 :   PetscCall(BVBiorthogonalizeGS(V,W,y));
     132         394 :   PetscCall(BVRestoreColumn(V,j,&y));
     133         394 :   PetscCall(BVGetColumn(W,j,&z));
     134         394 :   PetscCall(BVBiorthogonalizeGS(W,V,z));
     135         394 :   PetscCall(BVRestoreColumn(W,j,&z));
     136         394 :   V->k = ksavev;
     137         394 :   V->l = lsavev;
     138         394 :   W->k = ksavew;
     139         394 :   W->l = lsavew;
     140         394 :   PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
     141         394 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     142         394 :   PetscCall(PetscObjectStateIncrease((PetscObject)W));
     143         394 :   PetscFunctionReturn(PETSC_SUCCESS);
     144             : }
     145             : 
     146             : /*@
     147             :    BVBiorthonormalizeColumn - Bi-orthonormalize a column of two BV objects.
     148             : 
     149             :    Collective
     150             : 
     151             :    Input Parameters:
     152             : +  V - first basis vectors context
     153             : .  W - second basis vectors context
     154             : -  j - index of column to be bi-orthonormalized
     155             : 
     156             :    Output Parameters:
     157             : .  delta - (optional) value used for normalization
     158             : 
     159             :    Notes:
     160             :    This function first bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
     161             :    and V[0..j-1], respectively. Then, it scales the vectors with 1/delta, so
     162             :    that the resulting vectors satisfy W[j]'*V[j] = 1.
     163             : 
     164             :    Level: advanced
     165             : 
     166             : .seealso: BVOrthonormalizeColumn(), BVBiorthogonalizeColumn()
     167             : @*/
     168          73 : PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
     169             : {
     170          73 :   PetscScalar    alpha;
     171          73 :   PetscReal      deltat;
     172          73 :   PetscInt       ksavev,lsavev,ksavew,lsavew;
     173          73 :   Vec            y,z;
     174             : 
     175          73 :   PetscFunctionBegin;
     176          73 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     177          73 :   PetscValidHeaderSpecific(W,BV_CLASSID,2);
     178         219 :   PetscValidLogicalCollectiveInt(V,j,3);
     179          73 :   PetscValidType(V,1);
     180          73 :   BVCheckSizes(V,1);
     181          73 :   PetscValidType(W,2);
     182          73 :   BVCheckSizes(W,2);
     183          73 :   PetscCheckSameTypeAndComm(V,1,W,2);
     184          73 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     185          73 :   PetscCheck(j<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but V only has %" PetscInt_FMT " columns",j,V->m);
     186          73 :   PetscCheck(j<W->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but W only has %" PetscInt_FMT " columns",j,W->m);
     187          73 :   PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
     188          73 :   PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
     189          73 :   PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
     190          73 :   PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
     191             : 
     192             :   /* bi-orthogonalize */
     193          73 :   PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
     194          73 :   ksavev = V->k;
     195          73 :   lsavev = V->l;
     196          73 :   ksavew = W->k;
     197          73 :   lsavew = W->l;
     198          73 :   V->k = j;
     199          73 :   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
     200          73 :   W->k = j;
     201          73 :   W->l = -W->nc;
     202          73 :   PetscCall(BV_AllocateCoeffs(V));
     203          73 :   PetscCall(BV_AllocateCoeffs(W));
     204          73 :   PetscCall(BVGetColumn(V,j,&y));
     205          73 :   PetscCall(BVBiorthogonalizeGS(V,W,y));
     206          73 :   PetscCall(BVRestoreColumn(V,j,&y));
     207          73 :   PetscCall(BVGetColumn(W,j,&z));
     208          73 :   PetscCall(BVBiorthogonalizeGS(W,V,z));
     209          73 :   PetscCall(BVRestoreColumn(W,j,&z));
     210          73 :   V->k = ksavev;
     211          73 :   V->l = lsavev;
     212          73 :   W->k = ksavew;
     213          73 :   W->l = lsavew;
     214          73 :   PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
     215             : 
     216             :   /* scale */
     217          73 :   PetscCall(PetscLogEventBegin(BV_Scale,V,0,0,0));
     218          73 :   PetscCall(BVGetColumn(V,j,&y));
     219          73 :   PetscCall(BVGetColumn(W,j,&z));
     220          73 :   PetscCall(VecDot(z,y,&alpha));
     221          73 :   PetscCall(BVRestoreColumn(V,j,&y));
     222          73 :   PetscCall(BVRestoreColumn(W,j,&z));
     223          73 :   deltat = PetscSqrtReal(PetscAbsScalar(alpha));
     224          73 :   PetscUseTypeMethod(V,scale,j,1.0/PetscConj(alpha/deltat));
     225          73 :   PetscUseTypeMethod(W,scale,j,1.0/deltat);
     226          73 :   PetscCall(PetscLogEventEnd(BV_Scale,V,0,0,0));
     227          73 :   if (delta) *delta = deltat;
     228          73 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     229          73 :   PetscCall(PetscObjectStateIncrease((PetscObject)W));
     230          73 :   PetscFunctionReturn(PETSC_SUCCESS);
     231             : }