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 2339 : PetscErrorCode BVMatArnoldi(BV V,Mat A,Mat H,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
57 : {
58 2339 : PetscScalar *h;
59 2339 : const PetscScalar *a;
60 2339 : PetscInt j,ldh,rows,cols;
61 2339 : PetscBool lindep=PETSC_FALSE;
62 2339 : Vec buf;
63 :
64 2339 : PetscFunctionBegin;
65 2339 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
66 2339 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
67 7017 : PetscValidLogicalCollectiveInt(V,k,4);
68 2339 : PetscAssertPointer(m,5);
69 7017 : PetscValidLogicalCollectiveInt(V,*m,5);
70 2339 : PetscValidType(V,1);
71 2339 : BVCheckSizes(V,1);
72 2339 : PetscValidType(A,2);
73 2339 : PetscCheckSameComm(V,1,A,2);
74 :
75 2339 : 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 2339 : 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 2339 : PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
78 2339 : if (H) {
79 2339 : PetscValidHeaderSpecific(H,MAT_CLASSID,3);
80 2339 : PetscValidType(H,3);
81 2339 : PetscCheckTypeName(H,MATSEQDENSE);
82 2339 : PetscCall(MatGetSize(H,&rows,&cols));
83 2339 : PetscCall(MatDenseGetLDA(H,&ldh));
84 2339 : 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 2339 : 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 35925 : for (j=k;j<*m;j++) {
89 33615 : PetscCall(BVMatMultColumn(V,A,j));
90 33615 : 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 33595 : else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
92 33615 : if (PetscUnlikely(lindep)) {
93 29 : *m = j+1;
94 29 : break;
95 : }
96 : }
97 2339 : if (breakdown) *breakdown = lindep;
98 2339 : if (lindep) PetscCall(PetscInfo(V,"Arnoldi finished early at m=%" PetscInt_FMT "\n",*m));
99 :
100 2339 : if (H) {
101 2339 : PetscCall(MatDenseGetArray(H,&h));
102 2339 : PetscCall(BVGetBufferVec(V,&buf));
103 2339 : PetscCall(VecGetArrayRead(buf,&a));
104 33615 : for (j=k;j<*m-1;j++) PetscCall(PetscArraycpy(h+j*ldh,a+V->nc+(j+1)*(V->nc+V->m),j+2));
105 2339 : PetscCall(PetscArraycpy(h+(*m-1)*ldh,a+V->nc+(*m)*(V->nc+V->m),*m));
106 2339 : if (ldh>*m) h[(*m)+(*m-1)*ldh] = a[V->nc+(*m)+(*m)*(V->nc+V->m)];
107 2339 : PetscCall(VecRestoreArrayRead(buf,&a));
108 2339 : PetscCall(MatDenseRestoreArray(H,&h));
109 : }
110 :
111 2339 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
112 2339 : 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 1672 : PetscErrorCode BVMatLanczos(BV V,Mat A,Mat T,PetscInt k,PetscInt *m,PetscReal *beta,PetscBool *breakdown)
166 : {
167 1672 : PetscScalar *t;
168 1672 : const PetscScalar *a;
169 1672 : PetscReal *alpha,*betat;
170 1672 : PetscInt j,ldt,rows,cols,mincols=PetscDefined(USE_COMPLEX)?1:2;
171 1672 : PetscBool lindep=PETSC_FALSE;
172 1672 : Vec buf;
173 :
174 1672 : PetscFunctionBegin;
175 1672 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
176 1672 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
177 5016 : PetscValidLogicalCollectiveInt(V,k,4);
178 1672 : PetscAssertPointer(m,5);
179 5016 : PetscValidLogicalCollectiveInt(V,*m,5);
180 1672 : PetscValidType(V,1);
181 1672 : BVCheckSizes(V,1);
182 1672 : PetscValidType(A,2);
183 1672 : PetscCheckSameComm(V,1,A,2);
184 :
185 1672 : 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 1672 : 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 1672 : PetscCheck(*m>k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument m should be at least equal to k+1");
188 1672 : if (T) {
189 1672 : PetscValidHeaderSpecific(T,MAT_CLASSID,3);
190 1672 : PetscValidType(T,3);
191 1672 : PetscCheckTypeName(T,MATSEQDENSE);
192 1672 : PetscCall(MatGetSize(T,&rows,&cols));
193 1672 : PetscCall(MatDenseGetLDA(T,&ldt));
194 1672 : 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 1672 : 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 20434 : for (j=k;j<*m;j++) {
199 18782 : PetscCall(BVMatMultColumn(V,A,j));
200 18782 : 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 18762 : else PetscCall(BVOrthonormalizeColumn(V,j+1,PETSC_FALSE,beta,&lindep));
202 18782 : if (PetscUnlikely(lindep)) {
203 20 : *m = j+1;
204 20 : break;
205 : }
206 : }
207 1672 : if (breakdown) *breakdown = lindep;
208 1672 : if (lindep) PetscCall(PetscInfo(V,"Lanczos finished early at m=%" PetscInt_FMT "\n",*m));
209 :
210 1672 : if (T) {
211 1672 : PetscCall(MatDenseGetArray(T,&t));
212 1672 : alpha = (PetscReal*)t;
213 1672 : betat = alpha+ldt;
214 1672 : PetscCall(BVGetBufferVec(V,&buf));
215 1672 : PetscCall(VecGetArrayRead(buf,&a));
216 20454 : for (j=k;j<*m;j++) {
217 18782 : alpha[j] = PetscRealPart(a[V->nc+j+(j+1)*(V->nc+V->m)]);
218 18782 : betat[j] = PetscRealPart(a[V->nc+j+1+(j+1)*(V->nc+V->m)]);
219 : }
220 1672 : PetscCall(VecRestoreArrayRead(buf,&a));
221 1672 : PetscCall(MatDenseRestoreArray(T,&t));
222 : }
223 :
224 1672 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
225 1672 : PetscFunctionReturn(PETSC_SUCCESS);
226 : }
|
