LAM: SCALAPACK users?
Yoon Jae Ho
yoon at bh.kyungpook.ac.kr
Wed Aug 30 19:40:54 PDT 2000
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
More information about the Beowulf
mailing list