1 /*$Id: gmreig.c,v 1.28 2001/08/07 03:03:51 balay Exp $*/ 2 3 #include "src/ksp/ksp/impls/gmres/gmresp.h" 4 #include "petscblaslapack.h" 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "KSPComputeExtremeSingularValues_GMRES" 8 int KSPComputeExtremeSingularValues_GMRES(KSP ksp,PetscReal *emax,PetscReal *emin) 9 { 10 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 11 int n = gmres->it + 1,N = gmres->max_k + 2,ierr,lwork = 5*N,idummy = N,i; 12 PetscScalar *R = gmres->Rsvd,*work = R + N*N,sdummy; 13 PetscReal *realpart = gmres->Dsvd; 14 15 PetscFunctionBegin; 16 if (!n) { 17 *emax = *emin = 1.0; 18 PetscFunctionReturn(0); 19 } 20 /* copy R matrix to work space */ 21 ierr = PetscMemcpy(R,gmres->hh_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 22 23 /* zero below diagonal garbage */ 24 for (i=0; i<n; i++) { 25 R[i*N+i+1] = 0.0; 26 } 27 28 /* compute Singular Values */ 29 /* 30 The Cray math libraries on T3D/T3E, and Intel Math Kernel Libraries (MKL) for PCs do not 31 seem to have the DGESVD() lapack routines 32 */ 33 #if defined(PETSC_MISSING_LAPACK_GESVD) 34 SETERRQ(PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates."); 35 #else 36 #if !defined(PETSC_USE_COMPLEX) 37 LAgesvd_("N","N",&n,&n,R,&N,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&ierr); 38 #else 39 LAgesvd_("N","N",&n,&n,R,&N,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,realpart+N,&ierr); 40 #endif 41 if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in SVD Lapack routine"); 42 43 *emin = realpart[n-1]; 44 *emax = realpart[0]; 45 46 PetscFunctionReturn(0); 47 #endif 48 } 49 /* ------------------------------------------------------------------------ */ 50 /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */ 51 #if defined(PETSC_HAVE_ESSL) 52 #undef __FUNCT__ 53 #define __FUNCT__ "KSPComputeEigenvalues_GMRES" 54 int KSPComputeEigenvalues_GMRES(KSP ksp,int nmax,PetscReal *r,PetscReal *c,int *neig) 55 { 56 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 57 int n = gmres->it + 1,N = gmres->max_k + 1,ierr,lwork = 5*N; 58 int idummy = N,i,*perm,clen,zero; 59 PetscScalar *R = gmres->Rsvd; 60 PetscScalar *cwork = R + N*N,sdummy; 61 PetscReal *work,*realpart = gmres->Dsvd,*imagpart = realpart + N ; 62 63 PetscFunctionBegin; 64 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 65 *neig = n; 66 67 if (!n) { 68 PetscFunctionReturn(0); 69 } 70 /* copy R matrix to work space */ 71 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 72 73 /* compute eigenvalues */ 74 75 /* for ESSL version need really cwork of length N (complex), 2N 76 (real); already at least 5N of space has been allocated */ 77 78 ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr); 79 zero = 0; 80 LAgeev_(&zero,R,&N,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork); 81 ierr = PetscFree(work);CHKERRQ(ierr); 82 83 /* For now we stick with the convention of storing the real and imaginary 84 components of evalues separately. But is this what we really want? */ 85 ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr); 86 87 #if !defined(PETSC_USE_COMPLEX) 88 for (i=0; i<n; i++) { 89 realpart[i] = cwork[2*i]; 90 perm[i] = i; 91 } 92 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 93 for (i=0; i<n; i++) { 94 r[i] = cwork[2*perm[i]]; 95 c[i] = cwork[2*perm[i]+1]; 96 } 97 #else 98 for (i=0; i<n; i++) { 99 realpart[i] = PetscRealPart(cwork[i]); 100 perm[i] = i; 101 } 102 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 103 for (i=0; i<n; i++) { 104 r[i] = PetscRealPart(cwork[perm[i]]); 105 c[i] = PetscImaginaryPart(cwork[perm[i]]); 106 } 107 #endif 108 ierr = PetscFree(perm);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 #elif !defined(PETSC_USE_COMPLEX) 112 #undef __FUNCT__ 113 #define __FUNCT__ "KSPComputeEigenvalues_GMRES" 114 int KSPComputeEigenvalues_GMRES(KSP ksp,int nmax,PetscReal *r,PetscReal *c,int *neig) 115 { 116 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 117 int n = gmres->it + 1,N = gmres->max_k + 1,ierr,lwork = 5*N,idummy = N,i,*perm; 118 PetscScalar *R = gmres->Rsvd,*work = R + N*N; 119 PetscScalar *realpart = gmres->Dsvd,*imagpart = realpart + N,sdummy; 120 121 PetscFunctionBegin; 122 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 123 *neig = n; 124 125 if (!n) { 126 PetscFunctionReturn(0); 127 } 128 129 /* copy R matrix to work space */ 130 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 131 132 /* compute eigenvalues */ 133 #if defined(PETSC_MISSING_LAPACK_GEEV) 134 SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values."); 135 #else 136 LAgeev_("N","N",&n,R,&N,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&ierr); 137 #endif 138 if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in LAPACK routine"); 139 ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr); 140 for (i=0; i<n; i++) { perm[i] = i;} 141 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 142 for (i=0; i<n; i++) { 143 r[i] = realpart[perm[i]]; 144 c[i] = imagpart[perm[i]]; 145 } 146 ierr = PetscFree(perm);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 #else 150 #undef __FUNCT__ 151 #define __FUNCT__ "KSPComputeEigenvalues_GMRES" 152 int KSPComputeEigenvalues_GMRES(KSP ksp,int nmax,PetscReal *r,PetscReal *c,int *neig) 153 { 154 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 155 int n = gmres->it + 1,N = gmres->max_k + 1,ierr,lwork = 5*N,idummy = N,i,*perm; 156 PetscScalar *R = gmres->Rsvd,*work = R + N*N,*eigs = work + 5*N,sdummy; 157 158 PetscFunctionBegin; 159 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 160 *neig = n; 161 162 if (!n) { 163 PetscFunctionReturn(0); 164 } 165 /* copy R matrix to work space */ 166 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 167 168 /* compute eigenvalues */ 169 #if defined(PETSC_MISSING_LAPACK_GEEV) 170 SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values."); 171 #else 172 LAgeev_("N","N",&n,R,&N,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&ierr); 173 #endif 174 if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in LAPACK routine"); 175 ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr); 176 for (i=0; i<n; i++) { perm[i] = i;} 177 for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);} 178 ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr); 179 for (i=0; i<n; i++) { 180 r[i] = PetscRealPart(eigs[perm[i]]); 181 c[i] = PetscImaginaryPart(eigs[perm[i]]); 182 } 183 ierr = PetscFree(perm);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 #endif 187 188 189 190 191