LAM: SCALAPACK users?

[Camm Maguire]

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