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        1876 : static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
      45             : {
      46        1876 :   PetscFunctionBegin;
      47             :   /* h = W'*v */
      48        1876 :   PetscCall(BVDotVec(W,v,c));
      49             : 
      50             :   /* v = v - V h */
      51        1876 :   PetscCall(BVMultVec(V,-1.0,1.0,v,c));
      52             : 
      53        1876 :   PetscCall(BV_AddCoefficients(V,V->k,h,c));
      54        1876 :   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         994 : static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
      66             : {
      67         994 :   PetscScalar    *h,*c;
      68             : 
      69         994 :   PetscFunctionBegin;
      70         994 :   h = V->h;
      71         994 :   c = V->c;
      72         994 :   PetscCall(BV_CleanCoefficients(V,V->k,h));
      73        1932 :   PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
      74        1932 :   if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
      75         994 :   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         421 : PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
      97             : {
      98         421 :   PetscInt       ksavev,lsavev,ksavew,lsavew;
      99         421 :   Vec            y,z;
     100             : 
     101         421 :   PetscFunctionBegin;
     102         421 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     103         421 :   PetscValidHeaderSpecific(W,BV_CLASSID,2);
     104        1263 :   PetscValidLogicalCollectiveInt(V,j,3);
     105         421 :   PetscValidType(V,1);
     106         421 :   BVCheckSizes(V,1);
     107         421 :   PetscValidType(W,2);
     108         421 :   BVCheckSizes(W,2);
     109         421 :   PetscCheckSameTypeAndComm(V,1,W,2);
     110         421 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     111         421 :   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         421 :   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         421 :   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         421 :   PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
     115         421 :   PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
     116         421 :   PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
     117             : 
     118             :   /* bi-orthogonalize */
     119         421 :   PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
     120         421 :   ksavev = V->k;
     121         421 :   lsavev = V->l;
     122         421 :   ksavew = W->k;
     123         421 :   lsavew = W->l;
     124         421 :   V->k = j;
     125         421 :   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
     126         421 :   W->k = j;
     127         421 :   W->l = -W->nc;
     128         421 :   PetscCall(BV_AllocateCoeffs(V));
     129         421 :   PetscCall(BV_AllocateCoeffs(W));
     130         421 :   PetscCall(BVGetColumn(V,j,&y));
     131         421 :   PetscCall(BVBiorthogonalizeGS(V,W,y));
     132         421 :   PetscCall(BVRestoreColumn(V,j,&y));
     133         421 :   PetscCall(BVGetColumn(W,j,&z));
     134         421 :   PetscCall(BVBiorthogonalizeGS(W,V,z));
     135         421 :   PetscCall(BVRestoreColumn(W,j,&z));
     136         421 :   V->k = ksavev;
     137         421 :   V->l = lsavev;
     138         421 :   W->k = ksavew;
     139         421 :   W->l = lsavew;
     140         421 :   PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
     141         421 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     142         421 :   PetscCall(PetscObjectStateIncrease((PetscObject)W));
     143         421 :   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          76 : PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
     169             : {
     170          76 :   PetscScalar    alpha;
     171          76 :   PetscReal      deltat;
     172          76 :   PetscInt       ksavev,lsavev,ksavew,lsavew;
     173          76 :   Vec            y,z;
     174             : 
     175          76 :   PetscFunctionBegin;
     176          76 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     177          76 :   PetscValidHeaderSpecific(W,BV_CLASSID,2);
     178         228 :   PetscValidLogicalCollectiveInt(V,j,3);
     179          76 :   PetscValidType(V,1);
     180          76 :   BVCheckSizes(V,1);
     181          76 :   PetscValidType(W,2);
     182          76 :   BVCheckSizes(W,2);
     183          76 :   PetscCheckSameTypeAndComm(V,1,W,2);
     184          76 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     185          76 :   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          76 :   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          76 :   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          76 :   PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
     189          76 :   PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
     190          76 :   PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
     191             : 
     192             :   /* bi-orthogonalize */
     193          76 :   PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
     194          76 :   ksavev = V->k;
     195          76 :   lsavev = V->l;
     196          76 :   ksavew = W->k;
     197          76 :   lsavew = W->l;
     198          76 :   V->k = j;
     199          76 :   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
     200          76 :   W->k = j;
     201          76 :   W->l = -W->nc;
     202          76 :   PetscCall(BV_AllocateCoeffs(V));
     203          76 :   PetscCall(BV_AllocateCoeffs(W));
     204          76 :   PetscCall(BVGetColumn(V,j,&y));
     205          76 :   PetscCall(BVBiorthogonalizeGS(V,W,y));
     206          76 :   PetscCall(BVRestoreColumn(V,j,&y));
     207          76 :   PetscCall(BVGetColumn(W,j,&z));
     208          76 :   PetscCall(BVBiorthogonalizeGS(W,V,z));
     209          76 :   PetscCall(BVRestoreColumn(W,j,&z));
     210          76 :   V->k = ksavev;
     211          76 :   V->l = lsavev;
     212          76 :   W->k = ksavew;
     213          76 :   W->l = lsavew;
     214          76 :   PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
     215             : 
     216             :   /* scale */
     217          76 :   PetscCall(PetscLogEventBegin(BV_Scale,V,0,0,0));
     218          76 :   PetscCall(BVGetColumn(V,j,&y));
     219          76 :   PetscCall(BVGetColumn(W,j,&z));
     220          76 :   PetscCall(VecDot(z,y,&alpha));
     221          76 :   PetscCall(BVRestoreColumn(V,j,&y));
     222          76 :   PetscCall(BVRestoreColumn(W,j,&z));
     223          76 :   deltat = PetscSqrtReal(PetscAbsScalar(alpha));
     224          76 :   PetscUseTypeMethod(V,scale,j,1.0/PetscConj(alpha/deltat));
     225          76 :   PetscUseTypeMethod(W,scale,j,1.0/deltat);
     226          76 :   PetscCall(PetscLogEventEnd(BV_Scale,V,0,0,0));
     227          76 :   if (delta) *delta = deltat;
     228          76 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     229          76 :   PetscCall(PetscObjectStateIncrease((PetscObject)W));
     230          76 :   PetscFunctionReturn(PETSC_SUCCESS);
     231             : }