Example Scalapack pdsyev_

[murilo]

Hi all...

I need a simple example that call the subroutine in Scalapack pdsyev_ in C...
Somebody can help me??

Thanks..

Murilo

[bruce]

I post a similar question on how to call scalapack funtion in c. They gave me a good reply. Though that is a Lapack function but I think it might help. I am also trying to figure it out. My posted question is that how to get eigenvalues and eigenvectors. It is similar with using c to call fortran routine. Find some tutorials by search keywork "how to link c and fortran" in google. That might help.

Good luck

Bruce

[Julie]

Hi guys, here is a test code for pdsyev in C with my Makefile, Makefile.opts and include files.
test-pdsyev.c

Code: Select all

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/time.h>
#include "mpi.h"

#include "blas.h"
#include "blacs.h"
#include "scalapack.h"

extern void pdsyev_( char *jobz, char *uplo, int *n,
                double *a, int *ia, int *ja, int *desca, double *w, double *z, i                                                                                   nt *iz, int *jz, int *descz,
                double *work, int *lwork, int *info );


static int max( int a, int b ){
        if (a>b) return(a); else return(b);
}
static int min( int a, int b ){
        if (a<b) return(a); else return(b);
}

int main(int argc, char **argv) {
        int iam, nprocs;
        int myrank_mpi, nprocs_mpi;
        int ictxt, nprow, npcol, myrow, mycol;
        int np, nq, nb, n;
        int mpA, nqA;
        int i, j, k, info, itemp, seed, lwork, min_mn;
        int descA[9], descZ[9];
        double *A, *Z, *work, *W;
        int izero=0,ione=1;
        double mone=(-1.0e0),pone=(1.0e0),dzero=(0.0e0);
/**/
        double MPIt1, MPIt2, MPIelapsed, GFLOPS, GFLOPS_per_proc ;
        char jobz, uplo;
/**/
        MPI_Init( &argc, &argv);
        MPI_Comm_rank(MPI_COMM_WORLD, &myrank_mpi);
        MPI_Comm_size(MPI_COMM_WORLD, &nprocs_mpi);
/**/
        n = 100; nprow = 1; npcol = 1; nb = 64; jobz= 'V'; uplo='U';
        for( i = 1; i < argc; i++ ) {
                if( strcmp( argv[i], "-jobz" ) == 0 ) {
                        if (i+1<argc) {
                                if( strcmp( argv[i+1], "V" ) == 0 ){ jobz = 'V';                                                                                    i++; }
                                else if( strcmp( argv[i+1], "N" ) == 0 ){ jobz =                                                                                    'N'; i++; }
                                else if( strcmp( argv[i+1], "A" ) == 0 ){ jobz =                                                                                    'A'; i++; }
                                else printf(" ** warning: jobu should be set to                                                                                    V, N or A in the command line ** \n");
                        }
                        else
                                printf(" ** warning: jobu should be set to V, N                                                                                    or A in the command line ** \n");
                }
                if( strcmp( argv[i], "-n" ) == 0 ) {
                        n      = atoi(argv[i+1]);
                        i++;
                }
                if( strcmp( argv[i], "-p" ) == 0 ) {
                        nprow  = atoi(argv[i+1]);
                        i++;
                }
                if( strcmp( argv[i], "-q" ) == 0 ) {
                        npcol  = atoi(argv[i+1]);
                        i++;
                }
                if( strcmp( argv[i], "-nb" ) == 0 ) {
                        nb     = atoi(argv[i+1]);
                        i++;
                }
        }
/**/
        if (nb>n)
                nb = n;
        if (nprow*npcol>nprocs_mpi){
                if (myrank_mpi==0)
                        printf(" **** ERROR : we do not have enough processes av                                                                                   ailable to make a p-by-q process grid ***\n");
                        printf(" **** Bye-bye                                                                                                                                                            ***\n");
                MPI_Finalize(); exit(1);
        }
/**/
        Cblacs_pinfo( &iam, &nprocs ) ;
        Cblacs_get( -1, 0, &ictxt );
        Cblacs_gridinit( &ictxt, "Row", nprow, npcol );
        Cblacs_gridinfo( ictxt, &nprow, &npcol, &myrow, &mycol );
/**/
        if (iam==0)
                printf("\n");
/**/
/*
 *      if (iam==0)
 *              printf("\tn = %d\tnrhs = %d\t(%d,%d)\t%dx%d\n",n,nrhs,nprow,npco                                                                                   l,nb,nb);
 *      printf("Hello World, I am proc %d over %d for MPI, proc %d over %d for B                                                                                   LACS in position (%d,%d) in the process grid\n",
 *                      myrank_mpi,nprocs_mpi,iam,nprocs,myrow,mycol);
 */
/*
*
*     Work only the process in the process grid
*
*/
        if ((myrow < nprow)&(mycol < npcol)){
/*
*
*     Compute the size of the local matrices (thanks to numroc)
*
*/
                mpA    = numroc_( &n     , &nb, &myrow, &izero, &nprow );
                nqA    = numroc_( &n     , &nb, &mycol, &izero, &npcol );
/*
*
*     Allocate and fill the matrices A and B
*
*/
                A = (double *)calloc(mpA*nqA,sizeof(double)) ;
                if (A==NULL){ printf("error of memory allocation A on proc %dx%d                                                                                   \n",myrow,mycol); exit(0); }
                Z = (double *)calloc(mpA*nqA,sizeof(double)) ;
                if (Z==NULL){ printf("error of memory allocation VT on proc %dx%                                                                                   d\n",myrow,mycol); exit(0); }
                W = (double *)calloc(min_mn,sizeof(double)) ;
                if (W==NULL){ printf("error of memory allocation S on proc %dx%d                                                                                   \n",myrow,mycol); exit(0); }
/**/
                seed = iam*(mpA*nqA*2); srand(seed);
/**/
                k = 0;
                for (i = 0; i < mpA; i++) {
                        for (j = 0; j < nqA; j++) {
                                A[k] = ((double) rand()) / ((double) RAND_MAX) -                                                                                    0.5 ;
                                k++;
                        }
                }
/*
*
*     Initialize the array descriptor for the distributed matrices xA, U and VT
*
*/
                itemp = max( 1, mpA );
                descinit_( descA,  &n, &n, &nb, &nb, &izero, &izero, &ictxt, &it                                                                                   emp, &info );
                descinit_( descZ,  &n, &n, &nb, &nb, &izero, &izero, &ictxt, &it                                                                                   emp, &info );

                work = (double *)calloc(2,sizeof(double)) ;
                if (work==NULL){ printf("error of memory allocation for work on                                                                                    proc %dx%d (1st time)\n",myrow,mycol); exit(0); }
                lwork=-1;
                pdsyev_( &jobz, &uplo, &n, A, &ione, &ione, descA, W, Z, &ione,                                                                                    &ione, descZ, work, &lwork, &info );
                lwork= (int) work[0];
                free(work);
/**/
                work = (double *)calloc(lwork,sizeof(double)) ;
                if (work==NULL){ printf("error of memory allocation work on proc                                                                                    %dx%d\n",myrow,mycol); exit(0); }
/**/
                MPIt1 = MPI_Wtime();

                pdsyev_( &jobz, &uplo, &n, A, &ione, &ione, descA, W, Z, &ione,                                                                                    &ione, descZ, work, &lwork, &info );

                MPIt2 = MPI_Wtime();
                MPIelapsed=MPIt2-MPIt1;
/**/
                free(work);
                if ( iam==0 ){
                        printf("n=%d\t(%d,%d)\t%d\tjobz=%c\t%8.2fs \n",n,nprow,n                                                                                   pcol,nb,jobz,MPIelapsed);
                }
/**/
                free(W);
                free(Z);
                free(A);
        }
/*
*     Print ending messages
*/
        if ( iam==0 ){
                printf("\n");
        }
/**/
        Cblacs_gridexit( 0 );
        MPI_Finalize();
        exit(0);
}


Makefile

Code: Select all

include Makefile.opts

COMPILE = $(CC) $(CFLAGS) $(INCLUDEMPI) -c

%.o: %.c
        $(COMPILE) $*.c -o $@

SOURCES = test_psdyev.c
OBJECTS = $(SOURCES:%.c=%.o)

PDSYEV = test_pdsyev

all: $(PDSYEV)
pdsyev: $(PDSYEV)

LINK = $(F77) $(LDFLAGS)
LIBS = $(LIBSCALAPACK) $(LIBBLACS) $(LIBBLAS) $(LIBMPI)

test_pdsyev: test_pdsyev.o
        $(LINK) test_pdsyev.o $(LIBS) -o $@

clean:
        rm -f $(PDSYEV)

Makefile.opts

Code: Select all

CC = gcc
F77 = g77
CFLAGS = -O3 -Wall -I./include
INCLUDEMPI = -I/home/julie/opt/mpich-1.2.6/include
LDFLAGS = -O3
LIBSCALAPACK = /home/julie/lib/SCALAPACK_MPICH1/libscalapack.a
LIBBLACS = /home/julie/lib/BLACS_MPICH1/LIB/blacsF77init_MPI-LINUX-0.a /home/julie/lib/BLACS_MPICH1/LIB/blacs_MPI-LINUX-0.a /home/julie/lib/BLACS_MPICH1/LIB/blacsCinit_MPI-LINUX-0.a
/home/julie/lib/BLACS_MPICH1/LIB/blacsF77init_MPI-LINUX-0.a /home/julie/lib/BLACS_MPICH1/LIB/blacs_MPI-LINUX-0.a /home/julie/lib/BLACS_MPICH1/LIB/blacsCinit_MPI-LINUX-0.a
LIBBLAS = /usr/local/lib/libf77blas.a /usr/local/lib/libatlas.a
LIBMPI = /home/julie/opt/mpich-1.2.6/lib/libmpich.a

scalapack.h (that file is in a include directory)

Code: Select all

extern void   pdlawrite_( char **filenam, int *m, int *n, double *A, int *ia, int *ja, int *descA, int *irwrit, int *icwrit, double *work);
extern void   pdelset_( double *A, int *ia, int *ja, int *desca, double *alpha);
extern double pdlamch_( int *ictxt, char *cmach);
extern int    indxg2p_( int *indxglob, int *nb, int *iproc, int *isrcproc, int *nprocs);
extern int    indxg2l_( int *indxglob, int *nb, int *iproc, int *isrcproc, int *nprocs);
extern int    numroc_( int *n, int *nb, int *iproc, int *isrcproc, int *nprocs);
extern void   descinit_( int *desc, int *m, int *n, int *mb, int *nb, int *irsrc, int *icsrc,
                                int *ictxt, int *lld, int *info);
extern void   pdlaset_( char *uplo, int *m, int *n, double *alpha, double *beta, double *A, int *ia, int *ja, int *descA );
extern double pdlange_( char *norm, int *m, int *n, double *A, int *ia, int *ja, int *desca, double *work);
extern void   pdlacpy_( char *uplo, int *m, int *n, double *a, int *ia, int *ja, int *desca,
                                double *b, int *ib, int *jb, int *descb);
extern void   pdgesv_( int *n, int *nrhs, double *A, int *ia, int *ja, int *desca, int* ipiv,
                                double *B, int *ib, int *jb, int *descb, int *info);
extern void   pdgesvd_( char *jobu, char *jobvt, int *m, int *n, double *a, int *ia, int *ja, int *desca,
                                double *s, double *u, int *iu, int *ju, int *descu,
                                double *vt, int *ivt, int *jvt, int *descvt, double *work, int *lwork, int *info);
extern void   pdgemm_( char *TRANSA, char *TRANSB, int * M, int * N, int * K, double * ALPHA,
                                double * A, int * IA, int * JA, int * DESCA, double * B, int * IB, int * JB, int * DESCB,
                                double * BETA, double * C, int * IC, int * JC, int * DESCC );
extern int    indxg2p_( int *indxglob, int *nb, int *iproc, int *isrcproc, int *nprocs);

#ifdef F77_WITH_NO_UNDERSCORE
#define   numroc_      numroc
#define   descinit_    descinit
#define   pdlamch_     pdlamch
#define   pdlange_     pdlange
#define   pdlacpy_     pdlacpy
#define   pdgesv_      pdgesv
#define   pdgemm_      pdgemm
#define   indxg2p_     indxg2p
#endif

blacs.h (it is in the include directory)

Code: Select all

extern void   Cblacs_pinfo( int* mypnum, int* nprocs);
extern void   Cblacs_get( int context, int request, int* value);
extern int    Cblacs_gridinit( int* context, char * order, int np_row, int np_col);
extern void   Cblacs_gridinfo( int context, int*  np_row, int* np_col, int*  my_row, int*  my_col);
extern void   Cblacs_gridexit( int context);
extern void   Cblacs_exit( int error_code);


and last one:
blas.h

Code: Select all

extern void   dscal_( int *n, double *da, double *dx, int *incx);

Hope it helps
Julie

[bruce]

Hi, Julie:

Thank you for your code. It is very helpful. Based on your code, I have several questions.

1. How did you get the header files of blas.h, blacs.h and scalapack.h?
I find I have no these header files in the directories I installed. But I copy yours to compile your code. I am wondering if I need to write these header files by myself everytime and tailor them to some fortran functions? I copied your header files and tried to compiled it but I got slight error in the code.

2. In your Makefile.opts file, there is one line

CFLAGS = -O3 -Wall -I./include

I don't knw what this 'include' directory is. Which include directory it should be?

Thanks.

Bruce

[Julie]

Bruce,

Answer 1: Yes, I did the includes by hand. I am "old school".

Answer 2: the -I./include just means that my include files: blas.h, blacs.h, scalapack.h are in my include directory. With that, the compiler is able to find them.

If you have you include files where your source code is, just remove the -I./include, the compiler should look by himself in the current directory.

Julie

[itvasile]

Julie,
I have a question for you. I saw the code you posted and noticed that you have installed the scalapack/blacs in /home/julie. Why is that?
I'm new in this field, but I have some experience in Linux. I'm wondering because I saw in the scalapack install that it uses the $HOME variable for the current user to install the package. I don't like this kind of setting and I've installed it in /usr/local/scalapack (because I think that /home is for user documents only, not for installed programs that another user may want to use).
Thank you

[Julien Langou]

Well, ok, that's really not a big deal. Some prefer the installation to be in /usr/local, other wants the installation to be in a $HOME directory somewhere.

You are right if you are a sys-admin for sure you want your libscalapack.a to be in a location where users expect to find it and /usr/local/lib is certainly a good idea and /home/mysysadminlogin is certainly not a good idea. So once your installation is done, well, all you have to do is move the library file libscalpack.a from the $HOME directory to the /usr/local/lib directory. Not that hard. That's one option. It's kind of clean this way, in this sense that all the source files are not poluting the /usr directory.

The main reason for most of the people to like the default install to be in the $HOME directory is that they rarely have root permission and so they can not write in /usr .

In any case, I do not think that's a big issue. But keep us posted if you have a strong feeling in the matter.

Julien.

[itvasile]

Thanks for the answer. I'm indeed a sys-admin for those PCs and a programmer as well. So I have full control. My problem was that I saw too many posts with the installation in $HOME and I thought it was something like a "must" for scalapack (since I'm a newbie in this field). That's why I asked. Now I understood.

Thank you Julie for the makefiles. Yesterday I finished without any errors the linking and compiling of my parallel program. I didn't know there were that many dependancies.

John

[ncu571633]

Hi, Julien

There is a bug in your code.

W = (double *)calloc(min_mn,sizeof(double)) ;

min_mn is used before initialized.
It should be:
min_mn=min(mpA, mqA);

[Julien Langou]

A post from 2006 !

W should be of size N.

Moreover there is a typo in ScaLAPACK file:

Code: Select all

*  W       (global output) DOUBLE PRECISION array, dimension (N)
*          On normal exit, the first M entries contain the selected
*          eigenvalues in ascending order.


should be an N not an M.

Julien.