[Getdp] Q: Lanczos error
ohyeahq
ohyeahq at yahoo.co.jp
Thu Jul 1 09:29:46 CEST 2004
Christophe,
Thank you for your reply. I really like to use getdp for my problem (it's
TOO nice to be left unused!), so I tried to make lanczos work. It seems
ok,
although I haven't done extensive checking against Numerical Recipes
version. Here's my patch to Numeric/Lanczos.c and config.in. What I
did are:
(1) ran3() is replaced by GSL random number generator.
(2) GSL doesn't have hqr() analogue, so I use a LAPACK routine for
Hessenberg
eigensolver. Checking for LAPACK availability is added to config.in.
Relying on LAPACK is obviously minus in portability and performance (c to
fortran interfacing), but LAPACK routines are likely to be useful in future
extension of getdp (I hope).
Installing LAPACK is easy. For cygwin bash, I did something like this:
cd /usr/src
tar zxvf /path/to/lapack.tgz
cd LAPACK
cp INSTALL/make.inc.LINUX ./make.inc
make blaslib
make lapacklib
mv blas_LINUX.a /lib/libblas.a
mv lapack_LINUX.a /lib/liblapack.a
then
cd /path/to/getdp-src/
autoconf configure.in >configure
./configure
make
I tried a simple wave problem to compare with Numerical Recipes version.
They gave the same results, with minor difference such as different
eigenvalue ordering and eigenvector polarity. Performance seems similar
or slightly worse.
Hope this helps people who needs eigenvalues (like me).
Oh-Yeah? (Suzuki)
diff -c -r getdp-1.0.0_orig/Numeric/Lanczos.c getdp-1.0.0/Numeric/Lanczos.c
*** getdp-1.0.0_orig/Numeric/Lanczos.c Tue Jan 20 01:51:20 2004
--- getdp-1.0.0/Numeric/Lanczos.c Thu Jul 1 13:05:56 2004
***************
*** 177,187 ****
#if defined(HAVE_GSL)
! void Lanczos (struct DofData * DofData_P, int LanSize, List_T *LanSave,
double shift){
! Msg(ERROR, "Lanczos algorithm not available with the GSL");
}
! #else
/*
calcul des vecteurs propres d'une matrice de Hessenberg r馥lle
--- 177,309 ----
#if defined(HAVE_GSL)
! #include <gsl/gsl_rng.h>
!
! /* we use simple malloc for vectors & matrices (instead of
gsl_vector/matrix
! _set/get), so that convenient notation v[i] & a[i][j] can be used */
!
! /* base-1 index: numerical recipes convention */
! #define dvector(dummy, A) (double *)Malloc(sizeof(double) * ((A) + 1))
! #define dmatrix(dummy, A, dummy2, B) newmatrix((A) + 1, (B) + 1)
! #define free_dvector(A, dummy, dummy2) free(A)
! #define free_dmatrix(M, dummy, A, dummy2, dummy3) freematrix((M), (A))
!
! double **newmatrix(int nrow, int ncol) /* accessible as a[i][j] */
! {
! int i;
! double **a;
!
! a = (double **)Malloc(sizeof(double *) * nrow);
! if (a == NULL) {
! Msg(ERROR, "lanczos/newmatrix malloc error!");
! exit(1);
! }
! for (i = 0; i < nrow; i++) {
! a[i] = (double *)Malloc(sizeof(double) * ncol);
! if (a[i] == NULL) {
! while (--i >= 0) free(a[i]);
! free(a);
! Msg(ERROR, "lanczos/newmatrix malloc error!");
! exit(1);
! }
! }
! return a;
}
! void freematrix(double **a, int nrow)
! {
! int i;
! double **b;
!
! b = a;
! for (i = 0; i < nrow; i++) free(*b++);
! free(a);
! }
!
! double gsl_random(int isover) /* uniform random numbers in range [0.0,
1.0) */
! {
! double u;
! static const gsl_rng_type *T = NULL;
! static gsl_rng *r;
!
! if (T == NULL) { /* initialize */
! gsl_rng_env_setup();
! T = gsl_rng_default;
! r = gsl_rng_alloc(T);
! } else if (isover > 0) { /* finalize */
! gsl_rng_free(r);
! return 0.0;
! }
!
! u = gsl_rng_uniform(r);
! return u;
! }
!
! /* calc complex eigenvalues of Hessenberg form of real non-symmetrical
matrix:
! GSL doesn't have it, so we use LAPACK's DHSEQR routine */
!
! #if !defined(HAVE_LAPACK)
!
! void hqr(double **hin, int nhess, double *wrout, double *wiout)
! {
! Msg(ERROR, "Lanczos algorithm not available without LAPACK");
! }
!
! #else /* #if defined(HAVE_LAPACK) */
!
! #if defined(HAVE_UNDERSCORE)
! #define dhseqr_ dhseqr
! #endif
!
! void dhseqr_(char *JOB, char *COMPZ, int *N, int *ILO, int *IHI, double
*H,
! int *LDH, double WR[], double WI[], double Z[][], int *LDZ, double *WORK,
! int *LWORK, int *IRET);
!
! void hqr(double **hin, int nhess, double *wrout, double *wiout)
! {
! int n, ilo, ihi, ldh, ldz, lwork, info, i, j, k;
! double *h, *wr, *wi, *work, dummy[1][1];
! char job[] = "E", compz[] = "N";
!
! n = nhess;
! ilo = 1;
! ihi = n;
! ldh = n;
! ldz = n;
! lwork = n;
!
! wr = (double *)Malloc(sizeof(double) * n);
! wi = (double *)Malloc(sizeof(double) * n);
! work = (double *)Malloc(sizeof(double) * n);
!
! h = (double *)Malloc(sizeof(double) * (n * n));
!
! /* c to g77 interface: Fortran uses "column-major ordering" in arrays,
so
! transpose matrix and pack into h[] (or h[][]). This may or may not
! work with other systems. "cfortran.h" would be a better way */
! k = 0;
! for (i = 0; i < n; i++) {
! for (j = 0; j < n; j++) {
! h[k++] = hin[j + 1][i + 1]; /* hin: base-1 index */
! }
! }
!
! dhseqr_(job, compz, &n, &ilo, &ihi, h, &ldh, wr, wi, dummy, &ldz, work,
! &lwork, &info);
! if (info != 0) {
! Msg(BIGINFO, "lanczos/lapack dhseqr error = %d", info);
! }
!
! for (i = 0; i < n; i++) {
! wrout[i + 1] = wr[i]; /* base-1 */
! wiout[i + 1] = wi[i];
! }
! return;
! }
!
! #endif /* #if defined(HAVE_LAPACK) */
!
! #endif /* #if defined(HAVE_GSL) */
/*
calcul des vecteurs propres d'une matrice de Hessenberg r馥lle
***************
*** 258,264 ****
--- 380,388 ----
int i, j, k, ii, jj, NbrDof, Restart, ivmax, newsol = 0 ;
double dum, dum1, dum2;
double **Hes, **IMatrix, *diag, *wr, *wi, *err ;
+ #if !defined(HAVE_GSL)
long mun = -1 ;
+ #endif
struct Solution Solution_S ;
/* declaration pour les parametres de eigen.par */
***************
*** 384,391 ****
pourquoi ne pas essayer des 1 partout
Arnoldi-Tchebychev consiste ・am駘iorer le vecteur de d駱art
d'Arnoldi ・l'aide de Tchebychev */
for (i=0 ; i<NbrDof ; i++) LinAlg_SetDoubleInVector(ran3(&mun),
&Lan[0], i) ;
!
/* shifting */
if (fabs(shift) > 1.e-10)
LinAlg_AddMatrixProdMatrixDouble(K, M, -shift, K) ;
--- 508,520 ----
pourquoi ne pas essayer des 1 partout
Arnoldi-Tchebychev consiste ・am駘iorer le vecteur de d駱art
d'Arnoldi ・l'aide de Tchebychev */
+ #if defined(HAVE_GSL)
+ for (i = 0; i < NbrDof; i++) LinAlg_SetDoubleInVector(gsl_random(0),
+ &Lan[0], i);
+ gsl_random(1); /* finish it */
+ #else
for (i=0 ; i<NbrDof ; i++) LinAlg_SetDoubleInVector(ran3(&mun),
&Lan[0], i) ;
! #endif
/* shifting */
if (fabs(shift) > 1.e-10)
LinAlg_AddMatrixProdMatrixDouble(K, M, -shift, K) ;
***************
*** 797,800 ****
GetDP_End ;
}
- #endif /* #if !defined(HAVE_GSL) */
--- 926,928 ----
diff -c -r getdp-1.0.0_orig/configure.in getdp-1.0.0/configure.in
*** getdp-1.0.0_orig/configure.in Sun Mar 7 15:32:14 2004
--- getdp-1.0.0/configure.in Thu Jul 1 09:10:12 2004
***************
*** 221,226 ****
--- 221,234 ----
GETDP_LIBS="${GETDP_LIBS} -L${GSL_PREFIX}/lib -lgsl -lgslcblas"
C_FLAGS="${C_FLAGS} -I${GSL_PREFIX}/include"
fi
+ dnl Check for LAPACK (needed in Lanczos.c)
+ AC_CHECK_LIB(g2c,main)
+ AC_CHECK_LIB(blas,main)
+ AC_CHECK_LIB(lapack,main,LAPACK="yes",LAPACK="no")
+ if test "x${LAPACK}" = "xyes"; then
+ C_FLAGS="${C_FLAGS} -DHAVE_LAPACK"
+ GETDP_LIBS="${GETDP_LIBS} -llapack -lblas -lg2c"
+ fi
fi
fi
if test "x${GSL}" != "xyes"; then
__________________________________________________
Do You Yahoo!?
http://bb.yahoo.co.jp/