Actual source code: gmreig.c
2: #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
3: #include <petscblaslapack.h>
5: PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP ksp,PetscReal *emax,PetscReal *emin)
6: {
7: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
9: PetscInt n = gmres->it + 1,i,N = gmres->max_k + 2;
10: PetscBLASInt bn, bN,lwork, idummy,lierr;
11: PetscScalar *R = gmres->Rsvd,*work = R + N*N,sdummy = 0;
12: PetscReal *realpart = gmres->Dsvd;
15: PetscBLASIntCast(n,&bn);
16: PetscBLASIntCast(N,&bN);
17: PetscBLASIntCast(5*N,&lwork);
18: PetscBLASIntCast(N,&idummy);
19: if (n <= 0) {
20: *emax = *emin = 1.0;
21: return(0);
22: }
23: /* copy R matrix to work space */
24: PetscArraycpy(R,gmres->hh_origin,(gmres->max_k+2)*(gmres->max_k+1));
26: /* zero below diagonal garbage */
27: for (i=0; i<n; i++) R[i*N+i+1] = 0.0;
29: /* compute Singular Values */
30: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
31: #if !defined(PETSC_USE_COMPLEX)
32: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
33: #else
34: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,realpart+N,&lierr));
35: #endif
36: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr);
37: PetscFPTrapPop();
39: *emin = realpart[n-1];
40: *emax = realpart[0];
41: return(0);
42: }
44: PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
45: {
46: #if !defined(PETSC_USE_COMPLEX)
47: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
49: PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm;
50: PetscBLASInt bn, bN, lwork, idummy, l-1;
51: PetscScalar *R = gmres->Rsvd,*work = R + N*N;
52: PetscScalar *realpart = gmres->Dsvd,*imagpart = realpart + N,sdummy = 0;
55: PetscBLASIntCast(n,&bn);
56: PetscBLASIntCast(N,&bN);
57: PetscBLASIntCast(5*N,&lwork);
58: PetscBLASIntCast(N,&idummy);
59: if (nmax < n) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
60: *neig = n;
62: if (!n) return(0);
64: /* copy R matrix to work space */
65: PetscArraycpy(R,gmres->hes_origin,N*N);
67: /* compute eigenvalues */
68: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
69: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,R,&bN,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
70: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
71: PetscFPTrapPop();
72: PetscMalloc1(n,&perm);
73: for (i=0; i<n; i++) perm[i] = i;
74: PetscSortRealWithPermutation(n,realpart,perm);
75: for (i=0; i<n; i++) {
76: r[i] = realpart[perm[i]];
77: c[i] = imagpart[perm[i]];
78: }
79: PetscFree(perm);
80: #else
81: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
83: PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm;
84: PetscScalar *R = gmres->Rsvd,*work = R + N*N,*eigs = work + 5*N,sdummy;
85: PetscBLASInt bn,bN,lwork,idummy,l-1;
88: PetscBLASIntCast(n,&bn);
89: PetscBLASIntCast(N,&bN);
90: PetscBLASIntCast(5*N,&lwork);
91: PetscBLASIntCast(N,&idummy);
92: if (nmax < n) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
93: *neig = n;
95: if (!n) return(0);
97: /* copy R matrix to work space */
98: PetscArraycpy(R,gmres->hes_origin,N*N);
100: /* compute eigenvalues */
101: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
102: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,R,&bN,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&lierr));
103: if (lierr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine");
104: PetscFPTrapPop();
105: PetscMalloc1(n,&perm);
106: for (i=0; i<n; i++) perm[i] = i;
107: for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]);
108: PetscSortRealWithPermutation(n,r,perm);
109: for (i=0; i<n; i++) {
110: r[i] = PetscRealPart(eigs[perm[i]]);
111: c[i] = PetscImaginaryPart(eigs[perm[i]]);
112: }
113: PetscFree(perm);
114: #endif
115: return(0);
116: }
118: #if !defined(PETSC_USE_COMPLEX)
119: PetscErrorCode KSPComputeRitz_GMRES(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal *tetar,PetscReal *tetai)
120: {
121: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
123: PetscInt n = gmres->it + 1,N = gmres->max_k + 1,NbrRitz,nb=0;
124: PetscInt i,j,*perm;
125: PetscReal *H,*Q,*Ht; /* H Hessenberg Matrix and Q matrix of eigenvectors of H*/
126: PetscReal *wr,*wi,*modul; /* Real and imaginary part and modul of the Ritz values*/
127: PetscReal *SR,*work;
128: PetscBLASInt bn,bN,lwork,idummy;
129: PetscScalar *t,sdummy = 0;
132: /* n: size of the Hessenberg matrix */
133: if (gmres->fullcycle) n = N-1;
134: /* NbrRitz: number of (harmonic) Ritz pairs to extract */
135: NbrRitz = PetscMin(*nrit,n);
137: /* Definition of PetscBLASInt for lapack routines*/
138: PetscBLASIntCast(n,&bn);
139: PetscBLASIntCast(N,&bN);
140: PetscBLASIntCast(N,&idummy);
141: PetscBLASIntCast(5*N,&lwork);
142: /* Memory allocation */
143: PetscMalloc1(bN*bN,&H);
144: PetscMalloc1(bn*bn,&Q);
145: PetscMalloc1(lwork,&work);
146: PetscMalloc1(n,&wr);
147: PetscMalloc1(n,&wi);
149: /* copy H matrix to work space */
150: if (gmres->fullcycle) {
151: PetscArraycpy(H,gmres->hes_ritz,bN*bN);
152: } else {
153: PetscArraycpy(H,gmres->hes_origin,bN*bN);
154: }
156: /* Modify H to compute Harmonic Ritz pairs H = H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
157: if (!ritz) {
158: /* Transpose the Hessenberg matrix => Ht */
159: PetscMalloc1(bn*bn,&Ht);
160: for (i=0; i<bn; i++) {
161: for (j=0; j<bn; j++) {
162: Ht[i*bn+j] = H[j*bN+i];
163: }
164: }
165: /* Solve the system H^T*t = h^2_{m+1,m}e_m */
166: PetscCalloc1(bn,&t);
167: /* t = h^2_{m+1,m}e_m */
168: if (gmres->fullcycle) {
169: t[bn-1] = PetscSqr(gmres->hes_ritz[(bn-1)*bN+bn]);
170: } else {
171: t[bn-1] = PetscSqr(gmres->hes_origin[(bn-1)*bN+bn]);
172: }
173: /* Call the LAPACK routine dgesv to compute t = H^{-T}*t */
174: {
175: PetscBLASInt info;
176: PetscBLASInt nrhs = 1;
177: PetscBLASInt *ipiv;
178: PetscMalloc1(bn,&ipiv);
179: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn,&nrhs,Ht,&bn,ipiv,t,&bn,&info));
180: if (info) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Error while calling the Lapack routine DGESV");
181: PetscFree(ipiv);
182: PetscFree(Ht);
183: }
184: /* Now form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
185: for (i=0; i<bn; i++) H[(bn-1)*bn+i] += t[i];
186: PetscFree(t);
187: }
189: /* Compute (harmonic) Ritz pairs */
190: {
191: PetscBLASInt info;
192: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
193: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&bn,H,&bN,wr,wi,&sdummy,&idummy,Q,&bn,work,&lwork,&info));
194: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine");
195: }
196: /* sort the (harmonic) Ritz values */
197: PetscMalloc1(n,&modul);
198: PetscMalloc1(n,&perm);
199: for (i=0; i<n; i++) modul[i] = PetscSqrtReal(wr[i]*wr[i]+wi[i]*wi[i]);
200: for (i=0; i<n; i++) perm[i] = i;
201: PetscSortRealWithPermutation(n,modul,perm);
202: /* count the number of extracted Ritz or Harmonic Ritz pairs (with complex conjugates) */
203: if (small) {
204: while (nb < NbrRitz) {
205: if (!wi[perm[nb]]) nb += 1;
206: else nb += 2;
207: }
208: PetscMalloc1(nb*n,&SR);
209: for (i=0; i<nb; i++) {
210: tetar[i] = wr[perm[i]];
211: tetai[i] = wi[perm[i]];
212: PetscArraycpy(&SR[i*n],&(Q[perm[i]*bn]),n);
213: }
214: } else {
215: while (nb < NbrRitz) {
216: if (wi[perm[n-nb-1]] == 0) nb += 1;
217: else nb += 2;
218: }
219: PetscMalloc1(nb*n,&SR);
220: for (i=0; i<nb; i++) {
221: tetar[i] = wr[perm[n-nb+i]];
222: tetai[i] = wi[perm[n-nb+i]];
223: PetscArraycpy(&SR[i*n], &(Q[perm[n-nb+i]*bn]), n);
224: }
225: }
226: PetscFree(modul);
227: PetscFree(perm);
229: /* Form the Ritz or Harmonic Ritz vectors S=VV*Sr,
230: where the columns of VV correspond to the basis of the Krylov subspace */
231: if (gmres->fullcycle) {
232: for (j=0; j<nb; j++) {
233: VecZeroEntries(S[j]);
234: VecMAXPY(S[j],n,&SR[j*n],gmres->vecb);
235: }
236: } else {
237: for (j=0; j<nb; j++) {
238: VecZeroEntries(S[j]);
239: VecMAXPY(S[j],n,&SR[j*n],&VEC_VV(0));
240: }
241: }
242: *nrit = nb;
243: PetscFree(H);
244: PetscFree(Q);
245: PetscFree(SR);
246: PetscFree(wr);
247: PetscFree(wi);
248: return(0);
249: }
250: #endif