Sorry, forgot to attach the files....
--- Globe Trotter <itsme_410@xxxxxxxxx> wrote:
> Date: Sun, 6 Mar 2005 08:51:01 -0800 (PST)
> From: Globe Trotter <itsme_410@xxxxxxxxx>
> Subject: Serious Error: LAPACK and Fedora Core 3
> To: fedora <fedora-list@xxxxxxxxxx>
>
> 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 */
