Fwd: Serious Error: LAPACK and Fedora Core 3
Globe Trotter
itsme_410 at yahoo.com
Sun Mar 6 17:25:29 UTC 2005
OK,this seems a documented error on bugzilla. Fix seems an upgraded version of
lapack-3.0.28.
Though I may note that I used -g flag with the earlier 3.0.25 and still got the
hang....
As an aside, why not update the FC3 repository with the fix, instead of wasting
valuable time....I mean, I know, FC3 is bleeding-edge but once something is
documented and fixed, isn't it good to update with the fix and have users focus
on catching other bugs?
Thanks and best wishes!
--- Globe Trotter <itsme_410 at yahoo.com> wrote:
> Sorry, forgot to attach the files....
>
> --- Globe Trotter <itsme_410 at yahoo.com> wrote:
>
> > Date: Sun, 6 Mar 2005 08:51:01 -0800 (PST)
> > From: Globe Trotter <itsme_410 at yahoo.com>
> > Subject: Serious Error: LAPACK and Fedora Core 3
> > To: fedora <fedora-list at redhat.com>
> >
> > Dear all,
> >
> > I have been fighting with calling LAPACK from C and have come across an
> issue
> > which make eigenvalue calculations with LAPACK hang in Fedora Core 3. It
> may
> > be
> > noted that these go through fine in RHEL 3.
> >
> > I doubt this is an issue of calling LAPACK from C, but I will include the C
> > files and the header files.
> >
> > To compile,
> >
> > gcc demos.c mylapack.c -ansi -Wall -pedantic -llapack -lg2c
> >
> > Note that lapack needs to be installed. Also, the program compiles fine,
> and
> > calculates the inverse fine using LAPACK, on FC3, but hangs after entering
> > the
> > eigenvalue calculation. It goes to completion on RHEL 3.
> >
> > Can someone please suggest what can be done? Unless there is a strange bug
> in
> > FC3?
> >
> > By the way, the RHEL 3 box I tested on runs lapack-3.0.22, while the FC3
> box
> > runs lapack-3.0.25.
> >
> > Thanks and best wishes!
> >
> >
> >
> >
> >
> > __________________________________
> > Celebrate Yahoo!'s 10th Birthday!
> > Yahoo! Netrospective: 100 Moments of the Web
> > http://birthday.yahoo.com/netrospective/
> >
>
>
>
>
> __________________________________
> Celebrate Yahoo!'s 10th Birthday!
> Yahoo! Netrospective: 100 Moments of the Web
> http://birthday.yahoo.com/netrospective/> /*
> * Compile with:
> * cc try.c liblapack.a libblas.a -lg2c
> *
> * For a description of LINPACK's driver's routines, see:
> * http://www.netlib.org/lapack/lug/node25.html
> * */
>
> #include <stdio.h>
> #include <stdlib.h>
> #include "mylapack.h"
>
> static void print_matrix(const char *fmt, double *a, int m, int n)
> {
> int i, j;
>
> for (i=0; i<m; i++) {
> for (j=0; j<n; j++)
> printf(fmt, a[i+n*j]);
> putchar('\n');
> }
> }
>
> static void print_packed_matrix(const char *fmt, double *a, int n)
> {
> int i, j;
>
> for (i=0; i<n; i++) {
> for (j=0; j<=i; j++)
> printf(fmt, a[i+(j*(2*n-j-1))/2]);
> putchar('\n');
> }
> }
>
> static void print_vector(const char *fmt, double *b, int m)
> {
> int i;
>
> for (i=0; i<m; i++)
> printf(fmt, b[i]);
>
> putchar('\n');
> }
>
> /* LP_gen_solve() solves a general AX=B system */
> static void demo_gen_solve(int n)
> {
> double *a = NULL, *b = NULL;
> int *piv = NULL;
> int i, j;
> int info;
>
> a = malloc(n*n*sizeof *a);
> b = malloc(n*sizeof *b);
> piv = malloc(n*sizeof *piv);
>
> if (a==NULL || b==NULL || piv==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* make a sample non-symmetric matrix */
> for (j=0; j<n; j++)
> for (i=0; i<n; i++)
> a[i+n*j] = 1.0/(i+2*j+1);
>
> /* make a sample vector */
> for (i=0; i<n; i++)
> b[i] = 1.0;
>
> printf("A as a vector:\n");
> print_vector("%6.2f ", a, n*n);
>
> printf("A as a matrix:\n");
> print_matrix("%6.2f ", a, n, n);
>
> printf("B =\n");
> print_vector("%6.2f ", b, n);
>
> /* solve AX=B */
> info = LP_gen_solve(a, n, b, piv);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> printf("X =\n");
> print_vector("%6.2f ", b, n);
>
> cleanup:
> free(a);
> free(b);
> free(piv);
> }
>
> /* LP_sym_pos_solve() solves AX=B where A is symmetrix and positive-definite
> */
> static void demo_sym_pos_solve(int n)
> {
> double *a = NULL, *b = NULL;
> int i, j, k;
> int info;
>
> a = malloc(n*(n+1)/2*sizeof *a);
> b = malloc(n*sizeof *b);
>
> if (a==NULL || b==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* Make a sample symmetric matrix in the lower-triangular
> * packed format.
> * */
> for (j=0, k=0; j<n; j++)
> for (i=j; i<n; i++)
> a[k++] = 1.0/(i+j+1);
>
> /* make a sample vector */
> for (i=0; i<n; i++)
> b[i] = 1.0;
>
> printf("A =\n");
> print_packed_matrix("%6.2f ", a, n);
>
> printf("B =\n");
> print_vector("%6.2f ", b, n);
>
> /* solve AX=B */
> info = LP_sym_pos_solve(a, n, b);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> printf("X =\n");
> print_vector("%6.2f ", b, n);
>
> cleanup:
> free(a);
> free(b);
> }
>
> /* LP_sym_eigvals() finds the eigenvalues of a symmetrix matrix */
> static void demo_sym_eigvals(int n)
> {
> double *a = NULL, *w = NULL;
> int i, j;
> int info;
>
> a = malloc(n*n*sizeof *a);
> w = malloc(n*sizeof *w);
>
> if (a==NULL || w==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* Make a sample symmetric matrix */
> for (j=0; j<n; j++)
> for (i=0; i<n; i++)
> a[i+n*j] = 1.0/(i+j+1);
>
> printf("A =\n");
> print_matrix("%6.2f ", a, n, n);
>
> /* compute eigenvalues */
> info = LP_sym_eigvals(a, n, w);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> printf("eigenvalues =\n");
> print_vector("%g ", w, n);
>
> cleanup:
> free(a);
> free(w);
> }
>
> /* LP_sym_eigvecs() finds the eigenvalues and eigenvectors of
> * a symmetrix matrix
> * */
> static void demo_sym_eigvecs(int n)
> {
> double *a, *w, *z;
> int i, j;
> int info;
>
> a = malloc(n*n*sizeof *a);
> w = malloc(n*sizeof *w);
> z = malloc(n*n*sizeof *z);
>
> if (a==NULL || w==NULL || z==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* Make a sample symmetric matrix */
> for (j=0; j<n; j++)
> for (i=0; i<n; i++)
> a[i+n*j] = 1.0/(i+j+1);
>
> printf("A =\n");
> print_matrix("%6.2f ", a, n, n);
>
> /* compute eigenvalues and eigenvectors */
> info = LP_sym_eigvecs(a, n, w, z);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> printf("eigenvalues =\n");
> print_vector("%g ", w, n);
>
> printf("eigenvectors =\n");
> print_matrix("%12.8g ", z, n, n);
>
> cleanup:
> free(a);
> free(w);
> free(z);
> }
>
>
> int main(void)
> {
> int n = 3;
>
> printf("--- Solving AX=B using with LP_gen_solve() ---\n");
> demo_gen_solve(n);
>
> putchar('\n');
>
> printf("--- Solving AX=B with LP_sym_pos_solve()\n");
> demo_sym_pos_solve(n);
>
> putchar('\n');
>
> printf("--- Finding eigenvalues of A with LP_sym_eigvals() ---\n");
> demo_sym_eigvals(n);
>
> putchar('\n');
>
> printf("--- Finding eigenvalues and eigenvectors of A "
> "with LP_sym_eigvals() ---\n");
> demo_sym_eigvecs(n);
>
> return 0;
> }
> > #ifndef H_LAPACK_HEADERS_H
> #define H_LAPACK_HEADERS_H
>
> extern void dgesv_(int *n, int *nrhs, double *A, int *lda, int *piv,
> double *B, int *ldb, int *info);
>
> extern void dppsv_(char *uplo, int *n, int *nrhs, double *A,
> double *B, int *ldb, int *info);
>
> extern void dsyevr_(char *jobz, char *range, char *uplo,
> int *n, double *a, int *lda,
> double *vl, double *vu, int *il, int *iu,
> double *abstol, int *m, double *w,
> double *z, int *ldz, int *isuppz,
> double *work, int *lwork, int *iwork, int *liwork,
> int *info);
>
> #endif /* H_LAPACK_HEADERS_H */
> > #include <stdio.h>
> #include <stdlib.h>
> #include "lapack_headers.h"
> #include "mylapack.h"
>
> /* C front end to LAPACK's DGESV routine.
> *
> * Computes the solution of Ax=B, where A is an arbitrary real nxn matrix.
> *
> * The auxiliary vector piv returns the pivoting row exchange data.
> *
> * */
> int LP_gen_solve(double *a, int n, double *b, int *piv)
> {
> int nrhs = 1;
> int info;
>
> dgesv_(&n, &nrhs, a, &n, piv, b, &n, &info);
>
> return info;
> }
>
> /* C front end to LAPACK's DPPSV routine
> *
> * Computes the solution of AX=B where A is real symmetric and is stored
> * in the lower-triangular "packed" format.
> *
> */
> int LP_sym_pos_solve(double *a, int n, double *b)
> {
> int nrhs = 1;
> char uplo = 'L';
> int info;
>
> dppsv_(&uplo, &n, &nrhs, a, b, &n, &info);
>
> return info;
> }
>
> /* C front end to LINPACK's dsyevr_() routine.
> *
> * Calculates the eigenvalues of the nxn symmetric matrix A and returns
> * them in the vector w.
> *
> * */
> int LP_sym_eigvals(double *a, int n, double *w)
> {
> int lwork=-1, liwork=-1;
> double dx, *work = &dx; /* work points to a temporary cell */
> int ix, *iwork = &ix; /* iwork points to a temporary cell */
> double abstol = -1.0; /* force default */
> char jobz='N', range='A', uplo='L';
> int il, iu;
> double vl, vu;
> int info, m;
>
> /* Call dsyevr_() with lwork=-1 and liwork=-1 to query the
> * optimal sizes of the work and iwork arrays.
> * */
> dsyevr_(&jobz, &range, &uplo, &n, a, &n, &vl, &vu, &il, &iu,
> &abstol, &m, w, NULL, &n, NULL,
> work, &lwork, iwork, &liwork, &info);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> lwork = (int)*work;
> liwork = *iwork;
>
> /* allocate optimal sizes for work and iwork */
> work = malloc(lwork*sizeof *work);
> iwork = malloc(liwork*sizeof *iwork);
>
> if (work==NULL || iwork==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* now call dsyevr_() in earnest */
> dsyevr_(&jobz, &range, &uplo, &n, a, &n, &vl, &vu, &il, &iu,
> &abstol, &m, w, NULL, &n, NULL,
> work, &lwork, iwork, &liwork, &info);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> cleanup:
>
> free(work);
> free(iwork);
>
> return info;
> }
>
> /* C front end to LINPACK's dsyevr_() routine.
> *
> * Calculates the eigenvalues and eigenvectors of the nxn symmetric matrix A.
> * The eigenvalues are returned in the vector w.
> * The (orthonormal) eigenvectors are returned in the matrix z.
> * The ith column of z holds the eigenvector associated with w[i].
> *
> * */
> int LP_sym_eigvecs(double *a, int n, double *w, double *z)
> {
> int lwork=-1, liwork=-1;
> double dx, *work = &dx; /* work points to a temporary cell */
> int ix, *iwork = &ix; /* iwork points to a temporary cell */
> double abstol = -1.0; /* force default */
> int *isuppz;
> char jobz='V', range='A', uplo='L';
> int il, iu;
> double vl, vu;
> int info, m;
>
> /* Call dsyevr_() with lwork=-1 and liwork=-1 to query the
> * optimal sizes of the work and iwork arrays.
> * */
> dsyevr_(&jobz, &range, &uplo, &n, a, &n, &vl, &vu, &il, &iu,
> &abstol, &m, w, NULL, &n, NULL,
> work, &lwork, iwork, &liwork, &info);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> lwork = (int)*work;
> liwork = *iwork;
>
> /* allocate optimal sizes for work and iwork */
> work = malloc(lwork*sizeof *work);
> iwork = malloc(liwork*sizeof *iwork);
> isuppz = malloc(2*n*sizeof *isuppz);
>
> if (work==NULL || iwork==NULL || isuppz==NULL) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "out of memory!\n");
> goto cleanup;
> }
>
> /* now call dsyevr_() in earnest */
> dsyevr_(&jobz, &range, &uplo, &n, a, &n, &vl, &vu, &il, &iu,
> &abstol, &m, w, z, &n, isuppz,
> work, &lwork, iwork, &liwork, &info);
>
> if (info!=0) {
> fprintf(stderr, "trouble in file __FILE__, line __LINE__: "
> "info = %d\n", info);
> goto cleanup;
> }
>
> cleanup:
>
> free(work);
> free(iwork);
> free(isuppz);
>
> return info;
> }
>
> > #ifndef H_MYLAPACK_H
> #define H_MYLAPACK_H
>
> int LP_gen_solve(double *a, int n, double *b, int *piv);
> int LP_sym_pos_solve(double *a, int n, double *b);
> int LP_sym_eigvals(double *a, int n, double *w);
> int LP_sym_eigvecs(double *a, int n, double *w, double *z);
>
> #endif /* H_MYLAPACK_H */
> > --
> fedora-list mailing list
> fedora-list at redhat.com
> To unsubscribe: http://www.redhat.com/mailman/listinfo/fedora-list
__________________________________________________
Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around
http://mail.yahoo.com
More information about the fedora-list
mailing list