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 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr)); 32 #else 33 PetscStackCallBLAS("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 PetscCheckFalse(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 PetscStackCallBLAS("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 PetscStackCallBLAS("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 #if !defined(PETSC_USE_COMPLEX) 116 PetscErrorCode KSPComputeRitz_GMRES(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal *tetar,PetscReal *tetai) 117 { 118 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 119 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,NbrRitz,nb=0; 120 PetscInt i,j,*perm; 121 PetscReal *H,*Q,*Ht; /* H Hessenberg Matrix and Q matrix of eigenvectors of H*/ 122 PetscReal *wr,*wi,*modul; /* Real and imaginary part and modul of the Ritz values*/ 123 PetscReal *SR,*work; 124 PetscBLASInt bn,bN,lwork,idummy; 125 PetscScalar *t,sdummy = 0; 126 127 PetscFunctionBegin; 128 /* n: size of the Hessenberg matrix */ 129 if (gmres->fullcycle) n = N-1; 130 /* NbrRitz: number of (harmonic) Ritz pairs to extract */ 131 NbrRitz = PetscMin(*nrit,n); 132 133 /* Definition of PetscBLASInt for lapack routines*/ 134 PetscCall(PetscBLASIntCast(n,&bn)); 135 PetscCall(PetscBLASIntCast(N,&bN)); 136 PetscCall(PetscBLASIntCast(N,&idummy)); 137 PetscCall(PetscBLASIntCast(5*N,&lwork)); 138 /* Memory allocation */ 139 PetscCall(PetscMalloc1(bN*bN,&H)); 140 PetscCall(PetscMalloc1(bn*bn,&Q)); 141 PetscCall(PetscMalloc1(lwork,&work)); 142 PetscCall(PetscMalloc1(n,&wr)); 143 PetscCall(PetscMalloc1(n,&wi)); 144 145 /* copy H matrix to work space */ 146 if (gmres->fullcycle) { 147 PetscCall(PetscArraycpy(H,gmres->hes_ritz,bN*bN)); 148 } else { 149 PetscCall(PetscArraycpy(H,gmres->hes_origin,bN*bN)); 150 } 151 152 /* Modify H to compute Harmonic Ritz pairs H = H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */ 153 if (!ritz) { 154 /* Transpose the Hessenberg matrix => Ht */ 155 PetscCall(PetscMalloc1(bn*bn,&Ht)); 156 for (i=0; i<bn; i++) { 157 for (j=0; j<bn; j++) { 158 Ht[i*bn+j] = H[j*bN+i]; 159 } 160 } 161 /* Solve the system H^T*t = h^2_{m+1,m}e_m */ 162 PetscCall(PetscCalloc1(bn,&t)); 163 /* t = h^2_{m+1,m}e_m */ 164 if (gmres->fullcycle) { 165 t[bn-1] = PetscSqr(gmres->hes_ritz[(bn-1)*bN+bn]); 166 } else { 167 t[bn-1] = PetscSqr(gmres->hes_origin[(bn-1)*bN+bn]); 168 } 169 /* Call the LAPACK routine dgesv to compute t = H^{-T}*t */ 170 { 171 PetscBLASInt info; 172 PetscBLASInt nrhs = 1; 173 PetscBLASInt *ipiv; 174 PetscCall(PetscMalloc1(bn,&ipiv)); 175 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn,&nrhs,Ht,&bn,ipiv,t,&bn,&info)); 176 PetscCheck(!info,PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Error while calling the Lapack routine DGESV"); 177 PetscCall(PetscFree(ipiv)); 178 PetscCall(PetscFree(Ht)); 179 } 180 /* Now form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */ 181 for (i=0; i<bn; i++) H[(bn-1)*bn+i] += t[i]; 182 PetscCall(PetscFree(t)); 183 } 184 185 /* Compute (harmonic) Ritz pairs */ 186 { 187 PetscBLASInt info; 188 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 189 PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&bn,H,&bN,wr,wi,&sdummy,&idummy,Q,&bn,work,&lwork,&info)); 190 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine"); 191 } 192 /* sort the (harmonic) Ritz values */ 193 PetscCall(PetscMalloc1(n,&modul)); 194 PetscCall(PetscMalloc1(n,&perm)); 195 for (i=0; i<n; i++) modul[i] = PetscSqrtReal(wr[i]*wr[i]+wi[i]*wi[i]); 196 for (i=0; i<n; i++) perm[i] = i; 197 PetscCall(PetscSortRealWithPermutation(n,modul,perm)); 198 /* count the number of extracted Ritz or Harmonic Ritz pairs (with complex conjugates) */ 199 if (small) { 200 while (nb < NbrRitz) { 201 if (!wi[perm[nb]]) nb += 1; 202 else nb += 2; 203 } 204 PetscCall(PetscMalloc1(nb*n,&SR)); 205 for (i=0; i<nb; i++) { 206 tetar[i] = wr[perm[i]]; 207 tetai[i] = wi[perm[i]]; 208 PetscCall(PetscArraycpy(&SR[i*n],&(Q[perm[i]*bn]),n)); 209 } 210 } else { 211 while (nb < NbrRitz) { 212 if (wi[perm[n-nb-1]] == 0) nb += 1; 213 else nb += 2; 214 } 215 PetscCall(PetscMalloc1(nb*n,&SR)); 216 for (i=0; i<nb; i++) { 217 tetar[i] = wr[perm[n-nb+i]]; 218 tetai[i] = wi[perm[n-nb+i]]; 219 PetscCall(PetscArraycpy(&SR[i*n], &(Q[perm[n-nb+i]*bn]), n)); 220 } 221 } 222 PetscCall(PetscFree(modul)); 223 PetscCall(PetscFree(perm)); 224 225 /* Form the Ritz or Harmonic Ritz vectors S=VV*Sr, 226 where the columns of VV correspond to the basis of the Krylov subspace */ 227 if (gmres->fullcycle) { 228 for (j=0; j<nb; j++) { 229 PetscCall(VecZeroEntries(S[j])); 230 PetscCall(VecMAXPY(S[j],n,&SR[j*n],gmres->vecb)); 231 } 232 } else { 233 for (j=0; j<nb; j++) { 234 PetscCall(VecZeroEntries(S[j])); 235 PetscCall(VecMAXPY(S[j],n,&SR[j*n],&VEC_VV(0))); 236 } 237 } 238 *nrit = nb; 239 PetscCall(PetscFree(H)); 240 PetscCall(PetscFree(Q)); 241 PetscCall(PetscFree(SR)); 242 PetscCall(PetscFree(wr)); 243 PetscCall(PetscFree(wi)); 244 PetscFunctionReturn(0); 245 } 246 #endif 247