6
* $Id: cg.c,v 1.1.2.2 2007-07-05 15:29:56 vinod Exp $
19
# include <sys/types.h>
22
# include <sys/stat.h>
37
double *bvec,*dvec,*svec,*dmvec,*m_dvec,*amat,*xvec,*axvec,*rvec,*qvec;
40
int myfirstrow=0,mylastrow=0;
44
void read_and_create(int,char **);
45
void computeminverser(double *,double *, double *);
46
void computeminverse(double *,double *, int *,int *);
47
void finalize_arrays(int);
48
extern void acg_matvecmul(double *,double *,double *,int *,int *);
49
extern void acg_addvec(double *,double *,double *,double *, double *);
50
extern void acg_2addvec(double *,double *, double *,double *, double *,
51
double *, double *,double *, double *, double *, int *, int *);
52
extern double acg_ddot(double *,double *);
55
void conjugate_gradient(int nit,int dopreconditioning)
58
double d_one=1.0,d_negone=-1.0;
59
double delta0=0.0,deltaold=0.0,deltanew=0.0;
60
double alpha=0.0,negalpha,beta,dtransposeq;
62
acg_matvecmul(amat,xvec,axvec,ridx,cidx); /* compute Ax */
64
acg_printvec("axvec",axvec);
65
acg_printvec("bvec",bvec);
68
/* TODO original call had too many arguments */
69
/* acg_addvec(&d_one,bvec,&d_negone,axvec,rvec,ridx,cidx); */ /* r=b-Ax */
70
acg_addvec(&d_one,bvec,&d_negone,axvec,rvec); /* r=b-Ax */
73
acg_printvec("rvec",rvec);
74
if(me==0)for(i=0;i<nz;i++)printf("\n%d:col[%d]=%d",me,i,cidx[i]);
78
if(dopreconditioning){
79
computeminverse(dmvec,amat,ridx,cidx);
80
/*acg_printvec("dmvec",dmvec);*/
81
computeminverser(dmvec,rvec,dvec);
82
/*acg_printvec("dvec",dvec);*/
84
printf("\nDoing preconditioning!\n");
87
if(me==0)memcpy(dvec,rvec,na*sizeof(double));
90
deltanew = acg_ddot(rvec,dvec); /* deltanew = r.r_tranpose */
92
deltanew = acg_ddot(rvec,rvec); /* deltanew = r.r_tranpose */
94
delta0 = deltanew; /* delta0 = deltanew */
96
if(me==0)printf("\n\tdelta0 is %f\n",delta0);
98
for(i=0;i<nit && deltanew>(epsilon*epsilon*delta0);i++){
99
acg_matvecmul(amat,dvec,qvec,ridx,cidx); /* q = ad */
101
dtransposeq=acg_ddot(dvec,qvec); /* compute d_transpose.q */
103
alpha = deltanew/dtransposeq; /* deltanew/(d_transpose.q) */
107
acg_matvecmul(amat,xvec,axvec,ridx,cidx);
108
/* x = x+ alpha.d*/ /* r=b-Ax*/
109
acg_2addvec(&d_one,xvec,&alpha,dvec,xvec,&d_one,bvec,
110
&d_negone,axvec,rvec,ridx,cidx);
115
negalpha = 0.0-alpha;
116
/* x = x+ alpha.d*/ /* r=r-alpha.q*/
117
acg_2addvec(&d_one,xvec,&alpha,dvec,xvec,&d_one,rvec,
118
&negalpha,qvec,rvec,ridx,cidx);
121
if(dopreconditioning)
122
computeminverser(dmvec,rvec,svec);
124
deltaold = deltanew; /* deltaold = deltanew*/
126
if(dopreconditioning)
127
deltanew = acg_ddot(svec,rvec); /* deltanew = r_transpose.r*/
129
deltanew = acg_ddot(rvec,rvec); /* deltanew = r_transpose.r*/
131
beta = deltanew/deltaold; /* beta = deltanew/deltaold*/
133
if(dopreconditioning) {
134
/* too many aguments in function call */
135
/* acg_addvec(&d_one,svec,&beta,dvec,dvec,ridx,cidx); */
137
acg_addvec(&d_one,svec,&beta,dvec,dvec);
139
/* too many aguments in function call */
140
/* acg_addvec(&d_one,rvec,&beta,dvec,dvec,ridx,cidx); */
142
acg_addvec(&d_one,rvec,&beta,dvec,dvec);
146
acg_printvec("xvec",xvec);
148
/* acg_printvec("xvec",xvec); */
151
if(me==0)printf("\n\tIteration:%d\tBeta:%0.4f\tDelta:%f\n",
157
void initialize_arrays(int dpc)
162
dvec[i]=axvec[i]=rvec[i]=qvec[i]=0;
163
if(dpc){dmvec[i]=svec[i]=0;}
168
void finalize_arrays(int dpc)
191
int main(int argc, char **argv)
193
int dopreconditioning=1;
196
armci_msg_init(&argc,&argv);
197
nproc = armci_msg_nproc();
199
ARMCI_Init(); /* initialize ARMCI */
201
if(me==0)printf("\n CONJUGATE GRADIENT EXAMPLE\n");
204
printf(" CORRECT USAGE IS:");
205
printf("\n\n <launch commands> cg.x na nz file");
206
printf("\n\n where:");
207
printf("\n\tna is array dimention (only square arrays supported)");
208
printf("\n\tnz is number of non-zeros");
209
printf("\n\tfile is either the input file or the word random");
210
printf("\n\t use the word random if you to use random input");
211
printf("\n\t input should be in row compressed format");
212
printf("\n\t file should have matrix a followed by row, col & b (Ax=b)");
213
printf("\n\t if file also has na and nz, pass them as 0's and the");
214
printf("\n\t program will read them from the file");
215
printf("\n\nexample usages are:");
216
printf("\n\tmpirun -np 4 ./ga_cg.x 5000 80000 /home/me/myinput.dat");
218
printf("\n\tmpirun -np 4 ./ga_cg.x 5000 80000 random\n\n");
222
armci_msg_finalize();
226
read_and_create(argc,argv);
228
if(me==0)printf("\nWarmup and initialization run");
230
initialize_arrays(dopreconditioning);
231
conjugate_gradient(1,dopreconditioning);
234
if(me==0)printf("\n\nStarting Conjugate Gradient ....");
235
initialize_arrays(dopreconditioning);
237
conjugate_gradient(30000/*2*/,dopreconditioning);
240
acg_matvecmul(amat,xvec,axvec,ridx,cidx);
241
if(me==0)printf("\n%d:in %d iterations time to solution=%f-%f ax and b in cg_output.out\n",me,niter,(time1-time0),time_get);
242
acg_matvecmul(amat,xvec,axvec,ridx,cidx);
245
fd = fopen("cg_output.out", "w");
247
fprintf(fd,"\n%d:%s[%d]=%f %s[%d]=%f",me,"bvec",i,bvec[i],"axvec",i,axvec[i]);
252
finalize_arrays(dopreconditioning);
255
if(me==0)printf("Terminating ..\n");
257
armci_msg_finalize();