Actual source code: bvkrylov.c
slepc-3.21.1 2024-04-26
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: /*@C
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: }