LAM: SCALAPACK users?

[Yoon Jae Ho]

please refer http://www.netlib.org/scalapack/slug/node26.html

In those example1 , please do it like this using mpich ( not LAM )

g77 -O -o example1 example1.f /root/SCALAPACK/scalapack_LINUX.a   /root/SCALAPACK/pblas_LINUX.a  
         /root/SCALAPACK/tools_LINUX.a  /root/BLACS/LIB/blacsF77init_MPI-LINUX-0.a   /usr/local/blaslapack/blas.a
        /usr/local/mpich/build/LINUX/ch_p4/lib/libmpich.a  -lm  -lnsl

please refer the example1 from the above site .

from Yoon Jae Ho
http://ie.korea.ac.kr/~supercom/  Korea Beowulf
jhyoon at mail.posri.re.kr
yoon at bh.kyungpook.ac.kr
   

----- Original Message ----- 
From: Camm Maguire <camm at enhanced.com>
To: <william.s.yu at ieee.org>
Cc: <beowulf at beowulf.org>
Sent: Wednesday, August 30, 2000 11:42 PM
Subject: Re: LAM: SCALAPACK users?


> Greetings!  The following sample program (matrix definition ommitted),
> which implements an inverse iteration method to generate a particular
> eigenvector of a matrix, should clarify things, I hope.  Please write
> back if I can be of more assistance!
> 
> (P.S. One other thing I forgot to mention is that all fortran routines
> assume passed matrices are *column-major*.  The code below should
> clarify.)  
> 
> Take care,
> 
> =============================================================================
> main.c:
> =============================================================================
> #include <getopt.h>
> #include <math.h>
> #include <sys/time.h>
> #include <mpi.h>
> 
> #define TINY 1.0e-10
> 
> [snip]....
> 
> 
> static void
> mpiend(void) {
> 
>   int i;
> 
>   i=0;
>   blacs_exit__(&i);
>   
> 
> }
> 
> int
> main(int argc,char * argv[]) {
> 
>   int ic,nr=1,nc=1,mr,mc,N,nb=64,rsrc=0,mxlda,info,i,j,k,np,*iwork;
>   int desca[9],descaf[9],descb[9],descv[9],*ipiv,ione=1,lr,lc;
>   double **a,**af,**b,**v,*work,*x,*r,*c,rcond,ferr,berr,*q;
>   double xi,pp,pa,ps,pc,bx,t,l,l1;
> /*    struct timeval tv,tv1; */
>   int ch,debug=0,liwork=-1,lwork=-1,izero=0,itwo=2;
>   char *f,eq;
> 
>   
> 
>   MPI_Init(&argc,&argv);
>   MPI_Comm_size(MPI_COMM_WORLD,&np);
>   if (atexit(mpiend))
>     error("Can't setup mpiabort on exit\n");
> 
>   nc=nr=(int)sqrt(np);
> 
>   sl_init__(&ic,&nr,&nc);
>   blacs_gridinfo__(&ic,&nr,&nc,&mr,&mc);
> /*    Cblacs_gridinfo(ic,nr,nc,&mr,&mc); */
>   
>   i=0;
>   if (mr<0)
>     blacs_exit__(&i);
> 
>   pc=1.0;
>   
>   while ((ch=getopt(argc,argv,"x:i:p:a:s:c:n:d:"))!=EOF)
>     switch(ch) {
>     case 'x':
>       sscanf(optarg,"%lf",&bx);
>       break;
>     case 'i':
>       sscanf(optarg,"%lf",&xi);
>       break;
>     case 'p':
>       sscanf(optarg,"%lf",&pp);
>       break;
>     case 'a':
>       sscanf(optarg,"%lf",&pa);
>       break;
>     case 's':
>       sscanf(optarg,"%lf",&ps);
>       break;
>     case 'c':
>       sscanf(optarg,"%lf",&pc);
>       break;
>     case 'n':
>       sscanf(optarg,"%d",&N);
>       break;
>     case 'd':
>       sscanf(optarg,"%d",&debug);
>       break;
>     default:
>       break;
>     }
> 
>   mxlda=N/nc + N%nc;
>   if (mxlda%nb) mxlda+=(nb-mxlda%nb);
> 
>   lr=numroc_(&N,&nb,&mr,&rsrc,&nr);
>   lc=numroc_(&N,&nb,&mc,&rsrc,&nc);
> /*    lc=lr=mxlda; */
>   mem2(a,lc,lr);
>   mem2(af,lc,lr);
>   mem2(b,1,lr);
>   mem2(v,1,lr);
>   mem(x,N);
>   mem(ipiv,lr+nb);
>   mem(r,lr);
>   mem(c,lc);
>   mem(q,mxlda*mxlda);
> 
>   for (i=0;i<N;i++)
>     x[i]=bx+xi*i;
> 
>   l=1.0;
>   for (i=0;i<N;i++)
>     for (j=0;j<N;j++) {
>       int ni,nj;
> 
>       ni=i;
>       ni/=nb;
>       if (ni%nr!=mr)
> continue;
>       ni/=nr;
>       ni*=nb;
>       ni+=i%nb;
> 
>       nj=j;
>       nj/=nb;
>       if (nj%nc!=mc)
> continue;
>       nj/=nc;
>       nj*=nb;
>       nj+=j%nb;
> 
>       a[nj][ni]=kernel(x[i],x[j],0.0,pc*pp,pa/pc,ps/pc)*xi;
>       if (i==j)
> a[nj][ni]-=l;
> 
>     }
> 
>   descinit_(desca,&N,&N,&nb,&nb,&rsrc,&rsrc,&ic,&lr,&info);
>   descinit_(descaf,&N,&N,&nb,&nb,&rsrc,&rsrc,&ic,&lr,&info);
>   descinit_(descv,&N,&ione,&nb,&nb,&rsrc,&rsrc,&ic,&lr,&info);
>   descinit_(descb,&N,&ione,&nb,&nb,&rsrc,&rsrc,&ic,&lr,&info);
>   
>   if (!mc) {
>     randnor(b[0],lr,1.0);
>   
>     t=0.0;
>     for (i=0;i<lr;i++)
>       t+=b[0][i]*b[0][i];
>     dgsum2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>     dgebr2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>       
>     t=1.0/sqrt(t);
>     for (i=0;i<lr;i++)
>       b[0][i]*=t;
>   }
> 
> 
>   if (!mr && !mc) {
>     fprintf(stderr,"Starting ev search\n");
>     fflush(stderr);
>   }
> 
>   f="N";
>   k=-1;
>   pdgesvx_(f,"N",&N,&ione,a[0],&ione,&ione,desca,
>    af[0],&ione,&ione,descaf,ipiv,&eq,
>    r,c,b[0],&ione,&ione,descb,v[0],&ione,&ione,descv,
>    &rcond,&ferr,&berr,&t,&k,&j,&k,&info);
>   lwork=t;
>   mem(work,lwork);
>   liwork=j;
>   mem(iwork,liwork);
>   pdgesvx_(f,"N",&N,&ione,a[0],&ione,&ione,desca,
>    af[0],&ione,&ione,descaf,ipiv,&eq,
>    r,c,b[0],&ione,&ione,descb,v[0],&ione,&ione,descv,
>    &rcond,&ferr,&berr,work,&lwork,iwork,&liwork,&info);
>   if (debug && !mr && !mc)
>     fprintf(stderr,"Factor done: info %d rc %e fe %e be %e eq %c\n",
>     info,rcond,ferr,berr,eq);
>   if (info)
>     if (!mr && !mc)
>       error("Factorization failed: info %d rc %e fe %e be %e eq %c\n",
>     info,rcond,ferr,berr,eq);
>     else
>       return -1;
> 
>   f="F";
> 
>   for (i=0;i<200 && (!debug || i<debug);i++) {
> 
>     double t,t1[2];
>     int j;
> 
>     if (!mc) {
> 
>       for (j=0,t1[0]=t1[1]=0.0;j<lr;j++) {
> t1[0]+=v[0][j]*v[0][j];
> t1[1]+=v[0][j];
>       }
>       dgsum2d_(&ic,"C"," ",&itwo,&ione,t1,&itwo,&izero,&mc);
>       dgebr2d_(&ic,"C"," ",&itwo,&ione,t1,&itwo,&izero,&mc);
>       
>       t=1.0/sqrt(t1[0]);
>       t=t1[1]<0.0 ? -t : t;
>       
>       for (j=0;j<lr;j++)
> v[0][j]*=t;
> 
>       for (t=0.0,j=0;j<lr;j++)
> t+=(b[0][j]-v[0][j])*(b[0][j]-v[0][j]);
>       
>       dgsum2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>       dgebr2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>       
>       t=sqrt(t);
>     }
>     if (!mc)
>       memcpy(b[0],v[0],lr*sizeof(*b[0]));
> 
>     if (debug && !mc) {
> 
>       if (!mr)
> fprintf(stderr,"Result %d: diff %e\n",i,t);
> 
>       for (j=0;j<nr;j++) {
> if (j==mr) 
>   memcpy(q+j*mxlda,b[0],lr*sizeof(*b[0]));
> dgebr2d_(&ic,"C"," ",&mxlda,&ione,q+j*mxlda,&N,&j,&mc);
>       }
>       
>       if (!mr) {
> for (j=0;j<N;j++) {
>   int nj,nbb;
>   
>   nj=j;
>   nj/=nb;
>   nbb=nj%nr;
>   nj/=nr;
>   nj*=nb;
>   nj+=j%nb;
>   printf("%e %e\n",x[j],q[nbb*mxlda+nj]);
> }
> printf("\n\n");
>       }
>     }
> 
>     pdgesvx_(f,"N",&N,&ione,a[0],&ione,&ione,desca,
>      af[0],&ione,&ione,descaf,ipiv,&eq,
>      r,c,b[0],&ione,&ione,descb,v[0],&ione,&ione,descv,
>      &rcond,&ferr,&berr,work,&lwork,iwork,&liwork,&info);
>     if (debug && !mr && !mc)
>       fprintf(stderr,"Iteration %d: info %d rc %e fe %e be%e eq %c\n",
>       i,info,rcond,ferr,berr,eq);
> 
>     dgebr2d_(&ic,"R"," ",&ione,&ione,&t,&ione,&mr,&izero);
> 
>     if (t<TINY*N)
>       break;
> 
>   }
> 
>   if (!mc) {
> 
>     for (t=0.0,j=0;j<lr;j++)
>       t+=b[0][j]*v[0][j];
>     
>     dgsum2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>     dgebr2d_(&ic,"C"," ",&ione,&ione,&t,&ione,&izero,&mc);
>     
>     l1=l+1.0/t;
> 
>     if (!mr) {
>       fprintf(stderr,"done, %d iterations, eigenv %f\n",i,l1);
>       fflush(stderr);
>     }
> 
>     for (j=0;j<nr;j++) {
>       if (j==mr) 
> memcpy(q+j*mxlda,b[0],lr*sizeof(*b[0]));
>       dgebr2d_(&ic,"C"," ",&mxlda,&ione,q+j*mxlda,&N,&j,&mc);
>     }
>       
>     if (!mr)
>       for (i=0;i<N;i++) {
> int ni,nbb;
> 
> ni=i;
> ni/=nb;
> nbb=ni%nr;
> ni/=nr;
> ni*=nb;
> ni+=i%nb;
> printf("%e %e\n",x[i],q[nbb*mxlda+ni]);
>       }
>   }
> 
>   blacs_gridexit__(&ic);
> 
>   exit(0);
> 
> }
> =============================================================================
> Makefile:
> =============================================================================
> SC=-lscalapack-lam -lpblas-lam -ltools-lam -lredist-lam
> BL=-lblacsCinit-lam -lblacs-lam -lblacsCinit-lam
> BLAS=-lblas
> 
> LOADLIBES = $(SC) $(BL) $(BLAS) -lstd -lnum -lm -lmpi
> LINK=g77
> CFLAGS+=-I /usr/include/mpi
> 
> include $(HOME)/etc/Makefile
> 
> =============================================================================
> 
> "Horatio B. Bogbindero" <wyy at cersa.admu.edu.ph> writes:
> 
> > > 
> > > > are there any scalapack users in the list? i have a few questions to ask:
> > > > 
> > > 
> > > Here!
> > > 
> > > > -is there s C/C++ implementation of SCALAPACK? if not, can i call the
> > > > fortran SCALAPACK/PBLAS functions from C/C++?
> > > > 
> > > 
> > > No c interface, to use routines:
> > > 
> > > 1) add a '_' after routine name
> > > 2) Pass all arguments as pointers  (i.e. an arg of '1' would have to
> > > be passed as &i, with i having been set to 1.)
> > > 3) Link with g77
> > > 
> > hmmm. let me get this right. is it...
> > 
> > to call BLACS:
> > 
> > CALL BLACS_PINFO( IAM, NPROCS)
> > 
> > to c/c++ :
> > 
> > blacs_pinfo_(iam, nprocs);
> > 
> > to call a scalapack routine:
> > 
> > CALL PDGETRF( N( I ), N( I ), MEM( IPA ), 1, 1, DESCA, MEM( IPPIV ), INFO)
> > 
> > to c/c++:
> > 
> > pdgetrf_(n,n,mem1,1,1,desca,mem2,info);
> > 
> > ???
> > 
> > > > -are there any good scalapack documentation/manuals out there? the
> > > > scalapack site only feature some lawns but nothing like a users manual.
> > > > 
> > > 
> > > The faq and user's guide at netlib are both good.  They're included in
> > > the new Debian scalapack-doc package.
> > > 
> > how do i use the scalapack docs if they are a debian package and i use
> > redhat? is there a good printable version available like a pdf/ps version
> > of the scalapack-slug document at netlib.
> >  
> > ---------------------
> > william.s.yu at ieee.org
> >  
> > It's easier to fight for one's principles than to live up to them.
> >  
> > 
> > 
> > 
> 
> -- 
> Camm Maguire      camm at enhanced.com
> ==========================================================================
> "The earth is but one country, and mankind its citizens."  --  Baha'u'llah
> 
> _______________________________________________
> Beowulf mailing list
> Beowulf at beowulf.org
> http://www.beowulf.org/mailman/listinfo/beowulf