1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    BV routines related to Krylov decompositions
 12: */

 14: #include <slepc/private/bvimpl.h>

 16: /*@
 17:    BVMatArnoldi - Computes an Arnoldi factorization associated with a matrix.

 19:    Collective

 21:    Input Parameters:
 22: +  V - basis vectors context
 23: .  A - the matrix
 24: .  H - (optional) the upper Hessenberg matrix
 25: .  k - number of locked columns
 26: -  m - dimension of the Arnoldi basis, may be modified

 28:    Output Parameters:
 29: +  beta - (optional) norm of last vector before normalization
 30: -  breakdown - (optional) flag indicating that breakdown occurred

 32:    Notes:
 33:    Computes an m-step Arnoldi factorization for matrix A. The first k columns
 34:    are assumed to be locked and therefore they are not modified. On exit, the
 35:    following relation is satisfied

 37: $                    A * V - V * H = beta*v_m * e_m^T

 39:    where the columns of V are the Arnoldi vectors (which are orthonormal), H is
 40:    an upper Hessenberg matrix, e_m is the m-th vector of the canonical basis.
 41:    On exit, beta contains the norm of V[m] before normalization.

 43:    The breakdown flag indicates that orthogonalization failed, see
 44:    BVOrthonormalizeColumn(). In that case, on exit m contains the index of
 45:    the column that failed.

 47:    The values of k and m are not restricted to the active columns of V.

 49:    To create an Arnoldi factorization from scratch, set k=0 and make sure the
 50:    first column contains the normalized initial vector.

 52:    Level: advanced

 54: .seealso: BVMatLanczos(), BVSetActiveColumns(), BVOrthonormalizeColumn()
 55: @*/
 56: PetscErrorCode BVMatArnoldi(BV V,Mat A,Mat H,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
 57: {
 58:   PetscScalar       *h;
 59:   const PetscScalar *a;
 60:   PetscInt          j,ldh,rows,cols;
 61:   PetscBool         lindep=PETSC_FALSE;
 62:   Vec               buf;

 64:   PetscFunctionBegin;
 68:   PetscAssertPointer(m,5);
 71:   BVCheckSizes(V,1);
 73:   PetscCheckSameComm(V,1,A,2);

 75:   PetscCheck(k>=0 && k<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument k has wrong value %" PetscInt_FMT ", should be between 0 and %" PetscInt_FMT,k,V->m);
 76:   PetscCheck(*m>0 && *m<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m has wrong value %" PetscInt_FMT ", should be between 1 and %" PetscInt_FMT,*m,V->m);
 77:   PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
 78:   if (H) {
 81:     PetscCheckTypeName(H,MATSEQDENSE);
 82:     PetscCall(MatGetSize(H,&rows,&cols));
 83:     PetscCall(MatDenseGetLDA(H,&ldh));
 84:     PetscCheck(rows>=*m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Matrix H has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,rows,*m);
 85:     PetscCheck(cols>=*m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Matrix H has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,cols,*m);
 86:   }

 88:   for (j=k;j<*m;j++) {
 89:     PetscCall(BVMatMultColumn(V,A,j));
 90:     if (PetscUnlikely(j==V->N-1)) PetscCall(BV_OrthogonalizeColumn_Safe(V,j+1,NULL,beta,&lindep)); /* safeguard in case the full basis is requested */
 91:     else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
 92:     if (PetscUnlikely(lindep)) {
 93:       *m = j+1;
 94:       break;
 95:     }
 96:   }
 97:   if (breakdown) *breakdown = lindep;
 98:   if (lindep) PetscCall(PetscInfo(V,"Arnoldi finished early at m=%" PetscInt_FMT "\n",*m));

100:   if (H) {
101:     PetscCall(MatDenseGetArray(H,&h));
102:     PetscCall(BVGetBufferVec(V,&buf));
103:     PetscCall(VecGetArrayRead(buf,&a));
104:     for (j=k;j<*m-1;j++) PetscCall(PetscArraycpy(h+j*ldh,a+V->nc+(j+1)*(V->nc+V->m),j+2));
105:     PetscCall(PetscArraycpy(h+(*m-1)*ldh,a+V->nc+(*m)*(V->nc+V->m),*m));
106:     if (ldh>*m) h[(*m)+(*m-1)*ldh] = a[V->nc+(*m)+(*m)*(V->nc+V->m)];
107:     PetscCall(VecRestoreArrayRead(buf,&a));
108:     PetscCall(MatDenseRestoreArray(H,&h));
109:   }

111:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: /*@
116:    BVMatLanczos - Computes a Lanczos factorization associated with a matrix.

118:    Collective

120:    Input Parameters:
121: +  V - basis vectors context
122: .  A - the matrix
123: .  T - (optional) the tridiagonal matrix
124: .  k - number of locked columns
125: -  m - dimension of the Lanczos basis, may be modified

127:    Output Parameters:
128: +  beta - (optional) norm of last vector before normalization
129: -  breakdown - (optional) flag indicating that breakdown occurred

131:    Notes:
132:    Computes an m-step Lanczos factorization for matrix A, with full
133:    reorthogonalization. At each Lanczos step, the corresponding Lanczos
134:    vector is orthogonalized with respect to all previous Lanczos vectors.
135:    This is equivalent to computing an m-step Arnoldi factorization and
136:    exploting symmetry of the operator.

138:    The first k columns are assumed to be locked and therefore they are
139:    not modified. On exit, the following relation is satisfied

141: $                    A * V - V * T = beta*v_m * e_m^T

143:    where the columns of V are the Lanczos vectors (which are B-orthonormal),
144:    T is a real symmetric tridiagonal matrix, and e_m is the m-th vector of
145:    the canonical basis. On exit, beta contains the B-norm of V[m] before
146:    normalization. The T matrix is stored in a special way, its first column
147:    contains the diagonal elements, and its second column the off-diagonal
148:    ones. In complex scalars, the elements are stored as PetscReal and thus
149:    occupy only the first column of the Mat object. This is the same storage
150:    scheme used in matrix DS_MAT_T obtained with DSGetMat().

152:    The breakdown flag indicates that orthogonalization failed, see
153:    BVOrthonormalizeColumn(). In that case, on exit m contains the index of
154:    the column that failed.

156:    The values of k and m are not restricted to the active columns of V.

158:    To create a Lanczos factorization from scratch, set k=0 and make sure the
159:    first column contains the normalized initial vector.

161:    Level: advanced

163: .seealso: BVMatArnoldi(), BVSetActiveColumns(), BVOrthonormalizeColumn(), DSGetMat()
164: @*/
165: PetscErrorCode BVMatLanczos(BV V,Mat A,Mat T,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
166: {
167:   PetscScalar       *t;
168:   const PetscScalar *a;
169:   PetscReal         *alpha,*betat;
170:   PetscInt          j,ldt,rows,cols,mincols=PetscDefined(USE_COMPLEX)?1:2;
171:   PetscBool         lindep=PETSC_FALSE;
172:   Vec               buf;

174:   PetscFunctionBegin;
178:   PetscAssertPointer(m,5);
181:   BVCheckSizes(V,1);
183:   PetscCheckSameComm(V,1,A,2);

185:   PetscCheck(k>=0 || k<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument k has wrong value %" PetscInt_FMT ", should be between 0 and %" PetscInt_FMT,k,V->m);
186:   PetscCheck(*m>0 || *m<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m has wrong value %" PetscInt_FMT ", should be between 1 and %" PetscInt_FMT,*m,V->m);
187:   PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
188:   if (T) {
191:     PetscCheckTypeName(T,MATSEQDENSE);
192:     PetscCall(MatGetSize(T,&rows,&cols));
193:     PetscCall(MatDenseGetLDA(T,&ldt));
194:     PetscCheck(rows>=*m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Matrix T has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,rows,*m);
195:     PetscCheck(cols>=mincols,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Matrix T has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,cols,mincols);
196:   }

198:   for (j=k;j<*m;j++) {
199:     PetscCall(BVMatMultColumn(V,A,j));
200:     if (PetscUnlikely(j==V->N-1)) PetscCall(BV_OrthogonalizeColumn_Safe(V,j+1,NULL,beta,&lindep)); /* safeguard in case the full basis is requested */
201:     else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
202:     if (PetscUnlikely(lindep)) {
203:       *m = j+1;
204:       break;
205:     }
206:   }
207:   if (breakdown) *breakdown = lindep;
208:   if (lindep) PetscCall(PetscInfo(V,"Lanczos finished early at m=%" PetscInt_FMT "\n",*m));

210:   if (T) {
211:     PetscCall(MatDenseGetArray(T,&t));
212:     alpha = (PetscReal*)t;
213:     betat = alpha+ldt;
214:     PetscCall(BVGetBufferVec(V,&buf));
215:     PetscCall(VecGetArrayRead(buf,&a));
216:     for (j=k;j<*m;j++) {
217:       alpha[j] = PetscRealPart(a[V->nc+j+(j+1)*(V->nc+V->m)]);
218:       betat[j] = PetscRealPart(a[V->nc+j+1+(j+1)*(V->nc+V->m)]);
219:     }
220:     PetscCall(VecRestoreArrayRead(buf,&a));
221:     PetscCall(MatDenseRestoreArray(T,&t));
222:   }

224:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }