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