LAM: SCALAPACK users?
Camm Maguire
camm at enhanced.com
Wed Aug 30 07:42:56 PDT 2000
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
More information about the Beowulf
mailing list