mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-06 22:55:12 +00:00
Sample program now compiles and links correctly.
This commit is contained in:
Salvatore Filippone
+10
-11
@@ -1,41 +1,40 @@
|
||||
include ../../Make.inc
|
||||
include ../../../Make.inc
|
||||
#
|
||||
# Libraries used
|
||||
#
|
||||
LIBDIR=../../lib/
|
||||
INCLUDEDIR=../../include/
|
||||
LIBDIR=../../../lib/
|
||||
INCLUDEDIR=../../../include/
|
||||
HERE=.
|
||||
|
||||
FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBLAS_INCDIR)
|
||||
#PSBLAS_LIBS= -L$(PSBLAS_LIBDIR) -L$(LIBDIR) $(CPSBLAS_LIB) $(PSBLAS_LIB)
|
||||
# -lpsb_krylov_cbind -lpsb_prec_cbind -lpsb_base_cbind
|
||||
PSBC_LIBS= -L$(PSBLAS_CBIND_LIBDIR) \
|
||||
-lpsb_krylov_cbind -lpsb_prec_cbind -lpsb_base_cbind
|
||||
MLDC_LIBS=-L$(LIBDIR) -lmld_prec_cbind -lmld_prec -L$(MLDLIBDIR)
|
||||
PSBC_LIBS= -L$(PSBLAS_LIBDIR) -lpsb_cbind -lpsb_krylov -lpsb_prec
|
||||
MLDC_LIBS=-L$(LIBDIR) -lmld_cbind -lmld_prec
|
||||
#
|
||||
# Compilers and such
|
||||
#
|
||||
CCOPT= -g
|
||||
FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG).
|
||||
CINCLUDES=-I$(LIBDIR) -I$(PSBLAS_CBIND_INCDIR) -I$(INCLUDEDIR)
|
||||
CINCLUDES=-I$(LIBDIR) -I$(PSBLAS_INCDIR) -I$(INCLUDEDIR)
|
||||
|
||||
EXEDIR=./runs
|
||||
# 20110404 specifying UMFLIBS here is not portable!
|
||||
#UMFLIBS=-lumfpack -lamd -lcholmod -lcolamd -lcamd -lccolamd -L/usr/include/suitesparse
|
||||
#UMFFLAGS=-DHave_UMF_ -I/usr/include/suitesparse
|
||||
|
||||
all: mldec
|
||||
all: mldec
|
||||
|
||||
mldec: mldec.o
|
||||
$(MPF90) mldec.o -o mldec $(MLDC_LIBS) $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
|
||||
$(UMFLIBS) $(PSBLDLIBS) -lm -lgfortran
|
||||
$(MPFC) mldec.o -o mldec $(MLDC_LIBS) $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
|
||||
$(UMFLIBS) $(PSBLDLIBS) $(LDLIBS) -lm -lgfortran
|
||||
# \
|
||||
# -lifcore -lifcoremt -lguide -limf -lirc -lintlc -lcxaguard -L/opt/intel/fc/10.0.023/lib/ -lm
|
||||
|
||||
/bin/mv mldec $(EXEDIR)
|
||||
|
||||
.f90.o:
|
||||
$(MPF90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
|
||||
$(MPFC) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
|
||||
.c.o:
|
||||
$(MPCC) $(CCOPT) $(CINCLUDES) $(CDEFINES) -c $<
|
||||
|
||||
|
||||
+143
-139
@@ -76,151 +76,142 @@
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
|
||||
#include "psb_cbind.h"
|
||||
#include "psb_base_cbind.h"
|
||||
#include "mld_cbind.h"
|
||||
|
||||
|
||||
#define LINEBUFSIZE 1024
|
||||
#define NBMAX 20
|
||||
#define DUMPMATRIX 0
|
||||
|
||||
double a1(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(1.0/80.0);
|
||||
}
|
||||
double a2(double x, double y, double z)
|
||||
{
|
||||
return(20.0*y);
|
||||
return(1.0/80.0);
|
||||
}
|
||||
double a3(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(1.0/80.0);
|
||||
}
|
||||
double a4(double x, double y, double z)
|
||||
double c(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(0.0);
|
||||
}
|
||||
double b1(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(0.0/sqrt(3.0));
|
||||
}
|
||||
double b2(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(0.0/sqrt(3.0));
|
||||
}
|
||||
double b3(double x, double y, double z)
|
||||
{
|
||||
return(1.0);
|
||||
return(0.0/sqrt(3.0));
|
||||
}
|
||||
|
||||
int matgen(int ictxt, int ng,int idim,int vg[],psb_c_dspmat *ah,psb_c_descriptor *cdh,
|
||||
psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
|
||||
double g(double x, double y, double z)
|
||||
{
|
||||
int iam, np;
|
||||
int x, y, z, el,glob_row,i,info,ret;
|
||||
double gx, gy, gz, deltah;
|
||||
if (x == 1.0) {
|
||||
return(1.0);
|
||||
} else if (x == 0.0) {
|
||||
return( exp(-y*y-z*z));
|
||||
} else {
|
||||
return(0.0);
|
||||
}
|
||||
}
|
||||
|
||||
psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
|
||||
psb_c_dspmat *ah,psb_c_descriptor *cdh,
|
||||
psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
|
||||
{
|
||||
psb_i_t iam, np;
|
||||
psb_l_t ix, iy, iz, el,glob_row;
|
||||
psb_i_t i, k, info,ret;
|
||||
double x, y, z, deltah, sqdeltah, deltah2;
|
||||
double val[10*NBMAX], zt[NBMAX];
|
||||
int irow[10*NBMAX], icol[10*NBMAX];
|
||||
psb_l_t irow[10*NBMAX], icol[10*NBMAX];
|
||||
|
||||
info = 0;
|
||||
psb_c_info(ictxt,&iam,&np);
|
||||
deltah = (double) 1.0/(idim-1);
|
||||
|
||||
for (glob_row=1; glob_row<=ng; glob_row++) {
|
||||
|
||||
/* Check if I have to do something about this entry */
|
||||
if (vg[glob_row-1] == iam) {
|
||||
el=0;
|
||||
if ( (glob_row%(idim*idim)) == 0) {
|
||||
x = glob_row/(idim*idim);
|
||||
} else {
|
||||
x = glob_row/(idim*idim)+1;
|
||||
}
|
||||
if (((glob_row-(x-1)*idim*idim)%idim) == 0) {
|
||||
y = (glob_row-(x-1)*idim*idim)/idim;
|
||||
} else {
|
||||
y = (glob_row-(x-1)*idim*idim)/idim+1;
|
||||
}
|
||||
z = glob_row-(x-1)*idim*idim-(y-1)*idim;
|
||||
gx=x*deltah;
|
||||
gy=y*deltah;
|
||||
gz=z*deltah;
|
||||
zt[0] = 0.0;
|
||||
/* internal point: build discretization */
|
||||
/* term depending on (x-1,y,z) */
|
||||
|
||||
if (x==1) {
|
||||
val[el] = -b1(gx,gy,gz)-a1(gx,gy,gz);
|
||||
val[el] /= deltah*deltah;
|
||||
zt[0] = exp(-gy*gy-gz*gz)*(-val[el]);
|
||||
} else {
|
||||
val[el]=-b1(gx,gy,gz) -a1(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-2)*idim*idim+(y-1)*idim+(z);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y-1,z) */
|
||||
if (y==1) {
|
||||
val[el]=-b2(gx,gy,gz)-a2(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
zt[0] = exp(-gy*gy-gz*gz)*exp(-gx)*(-val[el]);
|
||||
} else {
|
||||
val[el]=-b2(gx,gy,gz)-a2(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-1)*idim*idim+(y-2)*idim+(z);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y,z-1)*/
|
||||
if (z==1) {
|
||||
val[el]=-b3(gx,gy,gz)-a3(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
zt[0] = exp(-gy*gy-gz*gz)*exp(-gx)*(-val[el]);
|
||||
} else {
|
||||
val[el]=-b3(gx,gy,gz)-a3(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-1)*idim*idim+(y-1)*idim+(z-1);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y,z)*/
|
||||
val[el]=2*b1(gx,gy,gz)+2*b2(gx,gy,gz)+2*b3(gx,gy,gz)
|
||||
+a1(gx,gy,gz)+a2(gx,gy,gz)+a3(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-1)*idim*idim+(y-1)*idim+(z);
|
||||
deltah = (double) 1.0/(idim+1);
|
||||
sqdeltah = deltah*deltah;
|
||||
deltah2 = 2.0* deltah;
|
||||
psb_c_set_index_base(0);
|
||||
for (i=0; i<nl; i++) {
|
||||
glob_row=vl[i];
|
||||
//if ((i%100000 == 0)||(i<10)) fprintf(stderr,"%d: generation loop at %d %ld \n",iam,i,glob_row);
|
||||
el=0;
|
||||
ix = glob_row/(idim*idim);
|
||||
iy = (glob_row-ix*idim*idim)/idim;
|
||||
iz = glob_row-ix*idim*idim-iy*idim;
|
||||
x=(ix+1)*deltah;
|
||||
y=(iy+1)*deltah;
|
||||
z=(iz+1)*deltah;
|
||||
zt[0] = 0.0;
|
||||
/* internal point: build discretization */
|
||||
/* term depending on (x-1,y,z) */
|
||||
val[el] = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2;
|
||||
if (ix==0) {
|
||||
zt[0] += g(0.0,y,z)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix-1)*idim*idim+(iy)*idim+(iz);
|
||||
el=el+1;
|
||||
/* term depending on (x,y,z+1) */
|
||||
if (z==idim) {
|
||||
val[el]=-b1(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
zt[0] = exp(-gy*gy-gz*gz)*exp(-gx)*(-val[el]);
|
||||
} else {
|
||||
val[el]=-b1(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-1)*idim*idim+(y-1)*idim+(z+1);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y+1,z) */
|
||||
if (y==idim) {
|
||||
val[el]=-b2(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
zt[0] = exp(-gy*gy-gz*gz)*exp(-gx)*(-val[el]);
|
||||
} else {
|
||||
val[el]=-b2(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x-1)*idim*idim+(y)*idim+(z);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x+1,y,z) */
|
||||
if (x<idim) {
|
||||
val[el]=-b3(gx,gy,gz);
|
||||
val[el] = val[el]/(deltah*deltah);
|
||||
icol[el]=(x)*idim*idim+(y-1)*idim+(z);
|
||||
el=el+1;
|
||||
}
|
||||
for (i=0; i<el; i++) irow[i]=glob_row;
|
||||
if ((ret=psb_c_dspins(el,irow,icol,val,ah,cdh))!=0)
|
||||
fprintf(stderr,"From psb_c_dspins: %d\n",ret);
|
||||
psb_c_dgeins(1,irow,zt,bh,cdh);
|
||||
zt[0]=0.0;
|
||||
psb_c_dgeins(1,irow,zt,xh,cdh);
|
||||
}
|
||||
/* term depending on (x,y-1,z) */
|
||||
val[el] = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2;
|
||||
if (iy==0) {
|
||||
zt[0] += g(x,0.0,z)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix)*idim*idim+(iy-1)*idim+(iz);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y,z-1)*/
|
||||
val[el]=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2;
|
||||
if (iz==0) {
|
||||
zt[0] += g(x,y,0.0)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix)*idim*idim+(iy)*idim+(iz-1);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y,z)*/
|
||||
val[el]=2.0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah + c(x,y,z);
|
||||
icol[el]=(ix)*idim*idim+(iy)*idim+(iz);
|
||||
el=el+1;
|
||||
/* term depending on (x,y,z+1) */
|
||||
val[el] = -a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2;
|
||||
if (iz==idim-1) {
|
||||
zt[0] += g(x,y,1.0)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix)*idim*idim+(iy)*idim+(iz+1);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x,y+1,z) */
|
||||
val[el] = -a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2;
|
||||
if (iy==idim-1) {
|
||||
zt[0] += g(x,1.0,z)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix)*idim*idim+(iy+1)*idim+(iz);
|
||||
el=el+1;
|
||||
}
|
||||
/* term depending on (x+1,y,z) */
|
||||
val[el] = -a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2;
|
||||
if (ix==idim-1) {
|
||||
zt[0] += g(1.0,y,z)*(-val[el]);
|
||||
} else {
|
||||
icol[el]=(ix+1)*idim*idim+(iy)*idim+(iz);
|
||||
el=el+1;
|
||||
}
|
||||
for (k=0; k<el; k++) irow[k]=glob_row;
|
||||
if ((ret=psb_c_dspins(el,irow,icol,val,ah,cdh))!=0)
|
||||
fprintf(stderr,"From psb_c_dspins: %d\n",ret);
|
||||
irow[0] = glob_row;
|
||||
psb_c_dgeins(1,irow,zt,bh,cdh);
|
||||
zt[0]=0.0;
|
||||
psb_c_dgeins(1,irow,zt,xh,cdh);
|
||||
}
|
||||
|
||||
if ((info=psb_c_cdasb(cdh))!=0) return(info);
|
||||
@@ -236,21 +227,22 @@ int matgen(int ictxt, int ng,int idim,int vg[],psb_c_dspmat *ah,psb_c_descriptor
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
int ictxt, iam, np;
|
||||
psb_i_t ictxt, iam, np;
|
||||
char methd[40], ptype[20], afmt[8], buffer[LINEBUFSIZE+1];
|
||||
int nparms;
|
||||
int idim,info,istop,itmax,itrace,irst,i,iter,ret;
|
||||
psb_i_t nparms;
|
||||
psb_i_t idim,info,istop,itmax,itrace,irst,iter,ret;
|
||||
mld_c_dprec *ph;
|
||||
psb_c_dspmat *ah;
|
||||
psb_c_dvector *bh, *xh, *rh;
|
||||
int *vg, ng, nb,nlr;
|
||||
psb_i_t nb,nlr, nl;
|
||||
psb_l_t i,ng, *vl, k;
|
||||
double t1,t2,eps,err;
|
||||
double *xv, *bv, *rv;
|
||||
double one=1.0, zero=0.0, res2;
|
||||
psb_c_SolverOptions options;
|
||||
psb_c_SolverOptions options;
|
||||
psb_c_descriptor *cdh;
|
||||
FILE *vectfile;
|
||||
|
||||
|
||||
ictxt = psb_c_init();
|
||||
psb_c_info(ictxt,&iam,&np);
|
||||
fprintf(stdout,"Initialization: am %d of %d\n",iam,np);
|
||||
@@ -291,18 +283,23 @@ int main(int argc, char *argv[])
|
||||
psb_c_barrier(ictxt);
|
||||
|
||||
cdh=psb_c_new_descriptor();
|
||||
psb_c_set_index_base(0);
|
||||
|
||||
/* Simple minded BLOCK data distribution */
|
||||
ng = idim*idim*idim;
|
||||
ng = ((psb_l_t) idim)*idim*idim;
|
||||
nb = (ng+np-1)/np;
|
||||
if ((vg=malloc(ng*sizeof(int)))==NULL) {
|
||||
nl = nb;
|
||||
if ( (ng -iam*nb) < nl) nl = ng -iam*nb;
|
||||
fprintf(stderr,"%d: Input data %d %ld %d %d\n",iam,idim,ng,nb, nl);
|
||||
if ((vl=malloc(nb*sizeof(psb_l_t)))==NULL) {
|
||||
fprintf(stderr,"On %d: malloc failure\n",iam);
|
||||
psb_c_abort(ictxt);
|
||||
}
|
||||
for (i=0; i<ng; i++) {
|
||||
vg[i] = i/nb;
|
||||
}
|
||||
if ((info=psb_c_cdall_vg(ng,vg,ictxt,cdh))!=0) {
|
||||
i = ((psb_l_t)iam) * nb;
|
||||
for (k=0; k<nl; k++)
|
||||
vl[k] = i+k;
|
||||
|
||||
if ((info=psb_c_cdall_vl(nl,vl,ictxt,cdh))!=0) {
|
||||
fprintf(stderr,"From cdall: %d\nBailing out\n",info);
|
||||
psb_c_abort(ictxt);
|
||||
}
|
||||
@@ -311,21 +308,26 @@ int main(int argc, char *argv[])
|
||||
xh = psb_c_new_dvector();
|
||||
rh = psb_c_new_dvector();
|
||||
ah = psb_c_new_dspmat();
|
||||
//fprintf(stderr,"From psb_c_new_dspmat: %p\n",ah);
|
||||
|
||||
/* Allocate mem space for sparse matrix and vectors */
|
||||
ret=psb_c_dspall(ah,cdh);
|
||||
//fprintf(stderr,"From psb_c_dspall: %d\n",ret);
|
||||
psb_c_dgeall(bh,cdh);
|
||||
psb_c_dgeall(xh,cdh);
|
||||
psb_c_dgeall(rh,cdh);
|
||||
|
||||
|
||||
|
||||
/* Matrix generation */
|
||||
if (matgen(ictxt,ng,idim,vg,ah,cdh,xh,bh,rh) != 0) {
|
||||
if (matgen(ictxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
|
||||
fprintf(stderr,"Error during matrix build loop\n");
|
||||
psb_c_abort(ictxt);
|
||||
}
|
||||
psb_c_barrier(ictxt);
|
||||
/* Set up the preconditioner */
|
||||
ph = mld_c_new_dprec();
|
||||
mld_c_dprecinit(ph,ptype,3);
|
||||
mld_c_dprecinit(ph,ptype);
|
||||
mld_c_dprecseti(ph,"SMOOTHER_SWEEPS",2);
|
||||
mld_c_dprecseti(ph,"SUB_FILLIN",1);
|
||||
mld_c_dprecsetc(ph,"COARSE_SOLVE","BJAC");
|
||||
@@ -341,15 +343,17 @@ int main(int argc, char *argv[])
|
||||
options.irst = irst;
|
||||
options.itrace = 1;
|
||||
options.istop = istop;
|
||||
|
||||
psb_c_seterraction_ret();
|
||||
t1=psb_c_wtime();
|
||||
ret=mld_c_dkrylov(methd,ah,ph,bh,xh,cdh,&options);
|
||||
t2=psb_c_wtime();
|
||||
iter = options.iter;
|
||||
err = options.err;
|
||||
fprintf(stderr,"From krylov: %d %lf, %d\n",iter,err,ret);
|
||||
|
||||
|
||||
//fprintf(stderr,"From krylov: %d %lf, %d %d\n",iter,err,ret,psb_c_get_errstatus());
|
||||
if (psb_c_get_errstatus() != 0) {
|
||||
psb_c_print_errmsg();
|
||||
}
|
||||
//fprintf(stderr,"After cleanup %d\n",psb_c_get_errstatus());
|
||||
/* Check 2-norm of residual on exit */
|
||||
psb_c_dgeaxpby(one,bh,zero,rh,cdh);
|
||||
psb_c_dspmm(-one,ah,xh,one,rh,cdh);
|
||||
@@ -363,11 +367,14 @@ int main(int argc, char *argv[])
|
||||
|
||||
}
|
||||
|
||||
#if 0
|
||||
bv = psb_c_dvect_get_cpy(bh);
|
||||
nlr=psb_c_cd_get_local_rows(cdh);
|
||||
#if DUMPATRIX
|
||||
psb_c_dmat_name_print(ah,"cbindmat.mtx");
|
||||
nlr = psb_c_cd_get_local_rows(cdh);
|
||||
bv = psb_c_dvect_get_cpy(bh);
|
||||
vectfile=fopen("cbindb.mtx","w");
|
||||
for (i=0;i<nlr; i++)
|
||||
fprintf(stdout,"RHS: %d %d %lf\n",iam,i,bv[i]);
|
||||
fprintf(vectfile,"%lf\n",bv[i]);
|
||||
fclose(vectfile);
|
||||
|
||||
|
||||
xv = psb_c_dvect_get_cpy(xh);
|
||||
@@ -401,7 +408,7 @@ int main(int argc, char *argv[])
|
||||
fprintf(stderr,"From cdfree: %d\nBailing out\n",info);
|
||||
psb_c_abort(ictxt);
|
||||
}
|
||||
fprintf(stderr,"pointer from cdfree: %p\n",cdh->descriptor);
|
||||
//fprintf(stderr,"pointer from cdfree: %p\n",cdh->descriptor);
|
||||
|
||||
/* Clean up object handles */
|
||||
free(ph);
|
||||
@@ -411,11 +418,8 @@ int main(int argc, char *argv[])
|
||||
free(cdh);
|
||||
|
||||
|
||||
fprintf(stderr,"program completed successfully\n");
|
||||
if (iam == 0) fprintf(stderr,"program completed successfully\n");
|
||||
|
||||
psb_c_barrier(ictxt);
|
||||
psb_c_exit(ictxt);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user