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 PetscErrorCode ierr; 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; 13 14 PetscFunctionBegin; 15 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 16 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 17 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 18 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 19 if (n <= 0) { 20 *emax = *emin = 1.0; 21 PetscFunctionReturn(0); 22 } 23 /* copy R matrix to work space */ 24 ierr = PetscArraycpy(R,gmres->hh_origin,(gmres->max_k+2)*(gmres->max_k+1));CHKERRQ(ierr); 25 26 /* zero below diagonal garbage */ 27 for (i=0; i<n; i++) R[i*N+i+1] = 0.0; 28 29 /* compute Singular Values */ 30 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 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 PetscCheckFalse(lierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr); 37 ierr = PetscFPTrapPop();CHKERRQ(ierr); 38 39 *emin = realpart[n-1]; 40 *emax = realpart[0]; 41 PetscFunctionReturn(0); 42 } 43 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; 48 PetscErrorCode ierr; 49 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm; 50 PetscBLASInt bn, bN, lwork, idummy, lierr = -1; 51 PetscScalar *R = gmres->Rsvd,*work = R + N*N; 52 PetscScalar *realpart = gmres->Dsvd,*imagpart = realpart + N,sdummy = 0; 53 54 PetscFunctionBegin; 55 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 56 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 57 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 58 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 59 PetscCheckFalse(nmax < n,PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 60 *neig = n; 61 62 if (!n) PetscFunctionReturn(0); 63 64 /* copy R matrix to work space */ 65 ierr = PetscArraycpy(R,gmres->hes_origin,N*N);CHKERRQ(ierr); 66 67 /* compute eigenvalues */ 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,R,&bN,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr)); 70 PetscCheckFalse(lierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr); 71 ierr = PetscFPTrapPop();CHKERRQ(ierr); 72 ierr = PetscMalloc1(n,&perm);CHKERRQ(ierr); 73 for (i=0; i<n; i++) perm[i] = i; 74 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 75 for (i=0; i<n; i++) { 76 r[i] = realpart[perm[i]]; 77 c[i] = imagpart[perm[i]]; 78 } 79 ierr = PetscFree(perm);CHKERRQ(ierr); 80 #else 81 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 82 PetscErrorCode ierr; 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,lierr = -1; 86 87 PetscFunctionBegin; 88 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 89 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 90 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 91 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 92 PetscCheckFalse(nmax < n,PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 93 *neig = n; 94 95 if (!n) PetscFunctionReturn(0); 96 97 /* copy R matrix to work space */ 98 ierr = PetscArraycpy(R,gmres->hes_origin,N*N);CHKERRQ(ierr); 99 100 /* compute eigenvalues */ 101 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 102 PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,R,&bN,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&lierr)); 103 PetscCheckFalse(lierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine"); 104 ierr = PetscFPTrapPop();CHKERRQ(ierr); 105 ierr = PetscMalloc1(n,&perm);CHKERRQ(ierr); 106 for (i=0; i<n; i++) perm[i] = i; 107 for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]); 108 ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr); 109 for (i=0; i<n; i++) { 110 r[i] = PetscRealPart(eigs[perm[i]]); 111 c[i] = PetscImaginaryPart(eigs[perm[i]]); 112 } 113 ierr = PetscFree(perm);CHKERRQ(ierr); 114 #endif 115 PetscFunctionReturn(0); 116 } 117 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; 122 PetscErrorCode ierr; 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; 130 131 PetscFunctionBegin; 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); 136 137 /* Definition of PetscBLASInt for lapack routines*/ 138 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 139 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 140 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 141 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 142 /* Memory allocation */ 143 ierr = PetscMalloc1(bN*bN,&H);CHKERRQ(ierr); 144 ierr = PetscMalloc1(bn*bn,&Q);CHKERRQ(ierr); 145 ierr = PetscMalloc1(lwork,&work);CHKERRQ(ierr); 146 ierr = PetscMalloc1(n,&wr);CHKERRQ(ierr); 147 ierr = PetscMalloc1(n,&wi);CHKERRQ(ierr); 148 149 /* copy H matrix to work space */ 150 if (gmres->fullcycle) { 151 ierr = PetscArraycpy(H,gmres->hes_ritz,bN*bN);CHKERRQ(ierr); 152 } else { 153 ierr = PetscArraycpy(H,gmres->hes_origin,bN*bN);CHKERRQ(ierr); 154 } 155 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 ierr = PetscMalloc1(bn*bn,&Ht);CHKERRQ(ierr); 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 ierr = PetscCalloc1(bn,&t);CHKERRQ(ierr); 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 ierr = PetscMalloc1(bn,&ipiv);CHKERRQ(ierr); 179 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn,&nrhs,Ht,&bn,ipiv,t,&bn,&info)); 180 PetscCheckFalse(info,PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Error while calling the Lapack routine DGESV"); 181 ierr = PetscFree(ipiv);CHKERRQ(ierr); 182 ierr = PetscFree(Ht);CHKERRQ(ierr); 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 ierr = PetscFree(t);CHKERRQ(ierr); 187 } 188 189 /* Compute (harmonic) Ritz pairs */ 190 { 191 PetscBLASInt info; 192 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 193 PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&bn,H,&bN,wr,wi,&sdummy,&idummy,Q,&bn,work,&lwork,&info)); 194 PetscCheckFalse(info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine"); 195 } 196 /* sort the (harmonic) Ritz values */ 197 ierr = PetscMalloc1(n,&modul);CHKERRQ(ierr); 198 ierr = PetscMalloc1(n,&perm);CHKERRQ(ierr); 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 ierr = PetscSortRealWithPermutation(n,modul,perm);CHKERRQ(ierr); 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 ierr = PetscMalloc1(nb*n,&SR);CHKERRQ(ierr); 209 for (i=0; i<nb; i++) { 210 tetar[i] = wr[perm[i]]; 211 tetai[i] = wi[perm[i]]; 212 ierr = PetscArraycpy(&SR[i*n],&(Q[perm[i]*bn]),n);CHKERRQ(ierr); 213 } 214 } else { 215 while (nb < NbrRitz) { 216 if (wi[perm[n-nb-1]] == 0) nb += 1; 217 else nb += 2; 218 } 219 ierr = PetscMalloc1(nb*n,&SR);CHKERRQ(ierr); 220 for (i=0; i<nb; i++) { 221 tetar[i] = wr[perm[n-nb+i]]; 222 tetai[i] = wi[perm[n-nb+i]]; 223 ierr = PetscArraycpy(&SR[i*n], &(Q[perm[n-nb+i]*bn]), n);CHKERRQ(ierr); 224 } 225 } 226 ierr = PetscFree(modul);CHKERRQ(ierr); 227 ierr = PetscFree(perm);CHKERRQ(ierr); 228 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 ierr = VecZeroEntries(S[j]);CHKERRQ(ierr); 234 ierr = VecMAXPY(S[j],n,&SR[j*n],gmres->vecb);CHKERRQ(ierr); 235 } 236 } else { 237 for (j=0; j<nb; j++) { 238 ierr = VecZeroEntries(S[j]);CHKERRQ(ierr); 239 ierr = VecMAXPY(S[j],n,&SR[j*n],&VEC_VV(0));CHKERRQ(ierr); 240 } 241 } 242 *nrit = nb; 243 ierr = PetscFree(H);CHKERRQ(ierr); 244 ierr = PetscFree(Q);CHKERRQ(ierr); 245 ierr = PetscFree(SR);CHKERRQ(ierr); 246 ierr = PetscFree(wr);CHKERRQ(ierr); 247 ierr = PetscFree(wi);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 #endif 251