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 routines related to Krylov decompositions
12 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 :
16 : /*@
17 : BVMatArnoldi - Computes an Arnoldi factorization associated with a matrix.
18 :
19 : Collective
20 :
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
27 :
28 : Output Parameters:
29 : + beta - (optional) norm of last vector before normalization
30 : - breakdown - (optional) flag indicating that breakdown occurred
31 :
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
36 :
37 : $ A * V - V * H = beta*v_m * e_m^T
38 :
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.
42 :
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.
46 :
47 : The values of k and m are not restricted to the active columns of V.
48 :
49 : To create an Arnoldi factorization from scratch, set k=0 and make sure the
50 : first column contains the normalized initial vector.
51 :
52 : Level: advanced
53 :
54 : .seealso: BVMatLanczos(), BVSetActiveColumns(), BVOrthonormalizeColumn()
55 : @*/
56 2269 : PetscErrorCode BVMatArnoldi(BV V,Mat A,Mat H,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
57 : {
58 2269 : PetscScalar *h;
59 2269 : const PetscScalar *a;
60 2269 : PetscInt j,ldh,rows,cols;
61 2269 : PetscBool lindep=PETSC_FALSE;
62 2269 : Vec buf;
63 :
64 2269 : PetscFunctionBegin;
65 2269 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
66 2269 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
67 6807 : PetscValidLogicalCollectiveInt(V,k,4);
68 2269 : PetscAssertPointer(m,5);
69 6807 : PetscValidLogicalCollectiveInt(V,*m,5);
70 2269 : PetscValidType(V,1);
71 2269 : BVCheckSizes(V,1);
72 2269 : PetscValidType(A,2);
73 2269 : PetscCheckSameComm(V,1,A,2);
74 :
75 2269 : 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 2269 : 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 2269 : PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
78 2269 : if (H) {
79 2269 : PetscValidHeaderSpecific(H,MAT_CLASSID,3);
80 2269 : PetscValidType(H,3);
81 2269 : PetscCheckTypeName(H,MATSEQDENSE);
82 2269 : PetscCall(MatGetSize(H,&rows,&cols));
83 2269 : PetscCall(MatDenseGetLDA(H,&ldh));
84 2269 : 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 2269 : 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 : }
87 :
88 34797 : for (j=k;j<*m;j++) {
89 32557 : PetscCall(BVMatMultColumn(V,A,j));
90 32557 : 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 32537 : else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
92 32557 : if (PetscUnlikely(lindep)) {
93 29 : *m = j+1;
94 29 : break;
95 : }
96 : }
97 2269 : if (breakdown) *breakdown = lindep;
98 2269 : if (lindep) PetscCall(PetscInfo(V,"Arnoldi finished early at m=%" PetscInt_FMT "\n",*m));
99 :
100 2269 : if (H) {
101 2269 : PetscCall(MatDenseGetArray(H,&h));
102 2269 : PetscCall(BVGetBufferVec(V,&buf));
103 2269 : PetscCall(VecGetArrayRead(buf,&a));
104 32557 : for (j=k;j<*m-1;j++) PetscCall(PetscArraycpy(h+j*ldh,a+V->nc+(j+1)*(V->nc+V->m),j+2));
105 2269 : PetscCall(PetscArraycpy(h+(*m-1)*ldh,a+V->nc+(*m)*(V->nc+V->m),*m));
106 2269 : if (ldh>*m) h[(*m)+(*m-1)*ldh] = a[V->nc+(*m)+(*m)*(V->nc+V->m)];
107 2269 : PetscCall(VecRestoreArrayRead(buf,&a));
108 2269 : PetscCall(MatDenseRestoreArray(H,&h));
109 : }
110 :
111 2269 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
112 2269 : PetscFunctionReturn(PETSC_SUCCESS);
113 : }
114 :
115 : /*@
116 : BVMatLanczos - Computes a Lanczos factorization associated with a matrix.
117 :
118 : Collective
119 :
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
126 :
127 : Output Parameters:
128 : + beta - (optional) norm of last vector before normalization
129 : - breakdown - (optional) flag indicating that breakdown occurred
130 :
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.
137 :
138 : The first k columns are assumed to be locked and therefore they are
139 : not modified. On exit, the following relation is satisfied
140 :
141 : $ A * V - V * T = beta*v_m * e_m^T
142 :
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().
151 :
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.
155 :
156 : The values of k and m are not restricted to the active columns of V.
157 :
158 : To create a Lanczos factorization from scratch, set k=0 and make sure the
159 : first column contains the normalized initial vector.
160 :
161 : Level: advanced
162 :
163 : .seealso: BVMatArnoldi(), BVSetActiveColumns(), BVOrthonormalizeColumn(), DSGetMat()
164 : @*/
165 2365 : PetscErrorCode BVMatLanczos(BV V,Mat A,Mat T,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
166 : {
167 2365 : PetscScalar *t;
168 2365 : const PetscScalar *a;
169 2365 : PetscReal *alpha,*betat;
170 2365 : PetscInt j,ldt,rows,cols,mincols=PetscDefined(USE_COMPLEX)?1:2;
171 2365 : PetscBool lindep=PETSC_FALSE;
172 2365 : Vec buf;
173 :
174 2365 : PetscFunctionBegin;
175 2365 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
176 2365 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
177 7095 : PetscValidLogicalCollectiveInt(V,k,4);
178 2365 : PetscAssertPointer(m,5);
179 7095 : PetscValidLogicalCollectiveInt(V,*m,5);
180 2365 : PetscValidType(V,1);
181 2365 : BVCheckSizes(V,1);
182 2365 : PetscValidType(A,2);
183 2365 : PetscCheckSameComm(V,1,A,2);
184 :
185 2365 : 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 2365 : 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 2365 : PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
188 2365 : if (T) {
189 2365 : PetscValidHeaderSpecific(T,MAT_CLASSID,3);
190 2365 : PetscValidType(T,3);
191 2365 : PetscCheckTypeName(T,MATSEQDENSE);
192 2365 : PetscCall(MatGetSize(T,&rows,&cols));
193 2365 : PetscCall(MatDenseGetLDA(T,&ldt));
194 2365 : 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 2365 : 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 : }
197 :
198 29032 : for (j=k;j<*m;j++) {
199 26687 : PetscCall(BVMatMultColumn(V,A,j));
200 26687 : 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 26667 : else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
202 26687 : if (PetscUnlikely(lindep)) {
203 20 : *m = j+1;
204 20 : break;
205 : }
206 : }
207 2365 : if (breakdown) *breakdown = lindep;
208 2365 : if (lindep) PetscCall(PetscInfo(V,"Lanczos finished early at m=%" PetscInt_FMT "\n",*m));
209 :
210 2365 : if (T) {
211 2365 : PetscCall(MatDenseGetArray(T,&t));
212 2365 : alpha = (PetscReal*)t;
213 2365 : betat = alpha+ldt;
214 2365 : PetscCall(BVGetBufferVec(V,&buf));
215 2365 : PetscCall(VecGetArrayRead(buf,&a));
216 29052 : for (j=k;j<*m;j++) {
217 26687 : alpha[j] = PetscRealPart(a[V->nc+j+(j+1)*(V->nc+V->m)]);
218 26687 : betat[j] = PetscRealPart(a[V->nc+j+1+(j+1)*(V->nc+V->m)]);
219 : }
220 2365 : PetscCall(VecRestoreArrayRead(buf,&a));
221 2365 : PetscCall(MatDenseRestoreArray(T,&t));
222 : }
223 :
224 2365 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
225 2365 : PetscFunctionReturn(PETSC_SUCCESS);
226 : }
|
