Actual source code: bvkrylov.c
slepc-3.18.1 2022-11-02
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 on V
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;
70: BVCheckSizes(V,1);
77: if (H) {
81: MatGetSize(H,&rows,&cols);
82: MatDenseGetLDA(H,&ldh);
85: }
87: for (j=k;j<*m;j++) {
88: BVMatMultColumn(V,A,j);
89: if (PetscUnlikely(j==V->N-1)) BV_OrthogonalizeColumn_Safe(V,j+1,NULL,beta,&lindep); /* safeguard in case the full basis is requested */
90: else BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep);
91: if (PetscUnlikely(lindep)) {
92: *m = j+1;
93: break;
94: }
95: }
96: if (breakdown) *breakdown = lindep;
97: if (lindep) PetscInfo(V,"Arnoldi finished early at m=%" PetscInt_FMT "\n",*m);
99: if (H) {
100: MatDenseGetArray(H,&h);
101: BVGetBufferVec(V,&buf);
102: VecGetArrayRead(buf,&a);
103: for (j=k;j<*m-1;j++) PetscArraycpy(h+j*ldh,a+V->nc+(j+1)*(V->nc+V->m),j+2);
104: PetscArraycpy(h+(*m-1)*ldh,a+V->nc+(*m)*(V->nc+V->m),*m);
105: if (ldh>*m) h[(*m)+(*m-1)*ldh] = a[V->nc+(*m)+(*m)*(V->nc+V->m)];
106: VecRestoreArrayRead(buf,&a);
107: MatDenseRestoreArray(H,&h);
108: }
110: PetscObjectStateIncrease((PetscObject)V);
111: return 0;
112: }
114: /*@C
115: BVMatLanczos - Computes a Lanczos factorization associated with a matrix.
117: Collective on V
119: Input Parameters:
120: + V - basis vectors context
121: . A - the matrix
122: . T - (optional) the tridiagonal matrix
123: . k - number of locked columns
124: - m - dimension of the Lanczos basis, may be modified
126: Output Parameters:
127: + beta - (optional) norm of last vector before normalization
128: - breakdown - (optional) flag indicating that breakdown occurred
130: Notes:
131: Computes an m-step Lanczos factorization for matrix A, with full
132: reorthogonalization. At each Lanczos step, the corresponding Lanczos
133: vector is orthogonalized with respect to all previous Lanczos vectors.
134: This is equivalent to computing an m-step Arnoldi factorization and
135: exploting symmetry of the operator.
137: The first k columns are assumed to be locked and therefore they are
138: not modified. On exit, the following relation is satisfied
140: $ A * V - V * T = beta*v_m * e_m^T
142: where the columns of V are the Lanczos vectors (which are B-orthonormal),
143: T is a real symmetric tridiagonal matrix, and e_m is the m-th vector of
144: the canonical basis. On exit, beta contains the B-norm of V[m] before
145: normalization. The T matrix is stored in a special way, its first column
146: contains the diagonal elements, and its second column the off-diagonal
147: ones. In complex scalars, the elements are stored as PetscReal and thus
148: occupy only the first column of the Mat object. This is the same storage
149: scheme used in matrix DS_MAT_T obtained with DSGetMat().
151: The breakdown flag indicates that orthogonalization failed, see
152: BVOrthonormalizeColumn(). In that case, on exit m contains the index of
153: the column that failed.
155: The values of k and m are not restricted to the active columns of V.
157: To create a Lanczos factorization from scratch, set k=0 and make sure the
158: first column contains the normalized initial vector.
160: Level: advanced
162: .seealso: BVMatArnoldi(), BVSetActiveColumns(), BVOrthonormalizeColumn(), DSGetMat()
163: @*/
164: PetscErrorCode BVMatLanczos(BV V,Mat A,Mat T,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
165: {
166: PetscScalar *t;
167: const PetscScalar *a;
168: PetscReal *alpha,*betat;
169: PetscInt j,ldt,rows,cols,mincols=PetscDefined(USE_COMPLEX)?1:2;
170: PetscBool lindep=PETSC_FALSE;
171: Vec buf;
179: BVCheckSizes(V,1);
186: if (T) {
190: MatGetSize(T,&rows,&cols);
191: MatDenseGetLDA(T,&ldt);
194: }
196: for (j=k;j<*m;j++) {
197: BVMatMultColumn(V,A,j);
198: if (PetscUnlikely(j==V->N-1)) BV_OrthogonalizeColumn_Safe(V,j+1,NULL,beta,&lindep); /* safeguard in case the full basis is requested */
199: else BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep);
200: if (PetscUnlikely(lindep)) {
201: *m = j+1;
202: break;
203: }
204: }
205: if (breakdown) *breakdown = lindep;
206: if (lindep) PetscInfo(V,"Lanczos finished early at m=%" PetscInt_FMT "\n",*m);
208: if (T) {
209: MatDenseGetArray(T,&t);
210: alpha = (PetscReal*)t;
211: betat = alpha+ldt;
212: BVGetBufferVec(V,&buf);
213: VecGetArrayRead(buf,&a);
214: for (j=k;j<*m;j++) {
215: alpha[j] = PetscRealPart(a[V->nc+j+(j+1)*(V->nc+V->m)]);
216: betat[j] = PetscRealPart(a[V->nc+j+1+(j+1)*(V->nc+V->m)]);
217: }
218: VecRestoreArrayRead(buf,&a);
219: MatDenseRestoreArray(T,&t);
220: }
222: PetscObjectStateIncrease((PetscObject)V);
223: return 0;
224: }