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: PetscFunctionReturn(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: . alpha - diagonal entries of tridiagonal matrix
123: . beta - subdiagonal entries of tridiagonal matrix
124: - k - number of locked columns
126: Input/Output Parameter:
127: . m - dimension of the Lanczos basis, may be modified
129: Output Parameter:
130: . breakdown - (optional) flag indicating that breakdown occurred
132: Notes:
133: Computes an m-step Lanczos factorization for matrix A, with full
134: reorthogonalization. At each Lanczos step, the corresponding Lanczos
135: vector is orthogonalized with respect to all previous Lanczos vectors.
136: This is equivalent to computing an m-step Arnoldi factorization and
137: exploting symmetry of the operator.
139: The first k columns are assumed to be locked and therefore they are
140: not modified. On exit, the following relation is satisfied
142: $ A * V - V * T = beta_m*v_m * e_m^T
144: where the columns of V are the Lanczos vectors (which are B-orthonormal),
145: T is a real symmetric tridiagonal matrix, and e_m is the m-th vector of
146: the canonical basis. The tridiagonal is stored as two arrays - alpha
147: contains the diagonal elements, beta the off-diagonal. On exit, the last
148: element of beta contains the B-norm of V[m] before normalization.
149: The basis V must have at least m+1 columns, while the arrays alpha and
150: beta must have space for at least m elements.
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()
164: @*/
165: PetscErrorCode BVMatLanczos(BV V,Mat A,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *m,PetscBool *breakdown)166: {
167: PetscScalar *a;
168: PetscInt j;
169: PetscBool lindep=PETSC_FALSE;
170: Vec buf;
180: BVCheckSizes(V,1);
188: for (j=k;j<*m;j++) {
189: BVMatMultColumn(V,A,j);
190: if (PetscUnlikely(j==V->N-1)) BV_OrthogonalizeColumn_Safe(V,j+1,NULL,beta+j,&lindep); /* safeguard in case the full basis is requested */
191: else BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta+j,&lindep);
192: if (PetscUnlikely(lindep)) {
193: *m = j+1;
194: break;
195: }
196: }
197: if (breakdown) *breakdown = lindep;
198: if (lindep) PetscInfo(V,"Lanczos finished early at m=%" PetscInt_FMT "\n",*m);
200: /* extract Hessenberg matrix from the BV buffer */
201: BVGetBufferVec(V,&buf);
202: VecGetArray(buf,&a);
203: for (j=k;j<*m;j++) alpha[j] = PetscRealPart(a[V->nc+j+(j+1)*(V->nc+V->m)]);
204: VecRestoreArray(buf,&a);
206: PetscObjectStateIncrease((PetscObject)V);
207: PetscFunctionReturn(0);
208: }