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