1 #if HAVE_CONFIG_H
2 # include "config.h"
3 #endif
4
5 #if HAVE_STDIO_H
6 # include <stdio.h>
7 #endif
8 #if HAVE_MATH_H
9 # include <math.h>
10 #endif
11 #if HAVE_STDLIB_H
12 # include <stdlib.h>
13 #endif
14 #if HAVE_SYS_TYPES_H
15 # include <sys/types.h>
16 #endif
17 #if HAVE_SYS_STAT_H
18 # include <sys/stat.h>
19 #endif
20 #if HAVE_FCNTL_H
21 # include <fcntl.h>
22 #endif
23 #if HAVE_STRING_H
24 # include <string.h>
25 #endif
26
27 #include "ga.h"
28 #include "macdecls.h"
29 #include "mp3.h"
30
31 extern int na;
32 extern int nz;
33 extern int dmvec,svec,bvec,dvec,m_dvec,amat,xvec,axvec,rvec,qvec,ridx,cidx;
34 extern int me, nproc;
35 extern int myfirstrow,mylastrow;
36 static int *columnmap,*allfirstrow,*alllastrow;
37 extern double *ga_vecptr;
38 extern int isvectormirrored;
39 static FILE *fd;
40
generate_random_file(int naa,int nnz)41 void generate_random_file(int naa,int nnz){
42 fd = fopen("randominput.dat", "w");
43 }
44
read_and_create(int argc,char ** argv)45 void read_and_create(int argc, char **argv)
46 {
47 int ri,i;
48 int ph;
49 double d_one=1.0;
50 double *a,*x;
51 int *icol, *irow;
52 #if 0
53 int dims[2];
54 #endif
55 int tmp1,idealelementsperproc;
56 int lo,hi;
57
58 na = atoi(argv[1]);
59 nz = atoi(argv[2]);
60
61 if(strncmp("random",argv[3],6)){
62 if(me==0){
63 fd = fopen(argv[3], "r");
64 if(fd==NULL)GA_Error("unable to open given file",0);
65 }
66 }
67 else{
68 if(na==0 || nz==0){
69 printf("\nERROR:exiting-no input file given and na or nz is 0");
70 fflush(stdout);
71 GA_Terminate();
72 MP_FINALIZE();
73 return;
74 }
75 if(me==0){
76 generate_random_file(na,nz);
77 fd = fopen("randominput.dat", "r");
78 }
79 }
80 if(me==0){
81 if(na==0)
82 fread(&na, sizeof(na), 1, fd);
83 if(nz==0)
84 fread(&nz, sizeof(nz), 1, fd);
85 printf("\nReading CG input\n");
86 printf("Number of rows: %d\n", na);
87 printf("Number of non-zeros: %d\n", nz);
88
89 a = (double *)malloc(sizeof(double)*nz);
90 icol = (int *)malloc(sizeof(int)*(nz+1));
91 x = (double *)malloc(sizeof(double)*(na+1));
92 irow = (int *)malloc(sizeof(int)*(na+1));
93
94 for (i = 0; i < na + 1; i++)
95 x[i] = 1.0;
96
97 fread(a, sizeof(double), nz, fd);
98 fread(irow, sizeof(int), na + 1, fd);
99 fread(icol, sizeof(int), nz + 1, fd);
100 fread(x, sizeof(double), na + 1, fd);
101
102 /* the c adjustment */
103 for (i = 0; i < na + 1; i++)
104 irow[i] -= 1;
105 for (i = 0; i < nz + 1; i++)
106 icol[i] -= 1;
107 /* MPI_Bcast(&nz,1,MPI_INT,0,MPI_COMM_WORLD); */
108 GA_Brdcst(&nz, sizeof(int), 0);
109 /* MPI_Bcast(&na,1,MPI_INT,0,MPI_COMM_WORLD); */
110 GA_Brdcst(&na, sizeof(int), 0);
111 }
112 else {
113 /* MPI_Bcast(&nz,1,MPI_INT,0,MPI_COMM_WORLD); */
114 GA_Brdcst(&nz, sizeof(int), 0);
115 /* MPI_Bcast(&na,1,MPI_INT,0,MPI_COMM_WORLD); */
116 GA_Brdcst(&na, sizeof(int), 0);
117 /*for now, others dont need to malloc really*/
118 a = (double *)malloc(sizeof(double)*nz);
119 icol = (int *)malloc(sizeof(int)*(nz+1));
120 x = (double *)malloc(sizeof(double)*(na+1));
121 irow = (int *)malloc(sizeof(int)*(na+1));
122 if(!a || !icol || !x || !irow)GA_Error("malloc failed in ga_cg",0);
123 }
124 allfirstrow = (int *)malloc(sizeof(int)*nproc);
125 alllastrow = (int *)malloc(sizeof(int)*nproc);
126 columnmap = (int *)malloc(sizeof(int)*nproc);
127 if(!allfirstrow || !alllastrow || !columnmap)GA_Error("malloc failed in ga_cg",0);
128
129 /*
130 * next decide who works on which rows, this will decide the
131 * distribution of a,d,r,q,x,and ax
132 */
133 /*create the mapping for all vectors, row matrix and column matrix*/
134 if(me==0){
135 idealelementsperproc = nz/nproc;
136 tmp1=0;
137 for(i=0;i<nproc;i++){
138 int elementsperproc=0;
139 allfirstrow[i]=tmp1;
140 for(ri=tmp1;ri<na;ri++,tmp1++){
141 elementsperproc+=(irow[ri+1]-irow[ri]);
142 if(elementsperproc>=idealelementsperproc){
143 if((elementsperproc-idealelementsperproc) >
144 idealelementsperproc-(elementsperproc-(irow[ri+1]-irow[ri]))){
145 alllastrow[i] = ri-1;
146 if((ri-1)<0)GA_Error("run on a smaller processor count",0);
147 tmp1--;
148 }
149 else{
150 alllastrow[i] = ri;
151 if(ri<0)GA_Error("run on a smaller processor count",0);
152 }
153 elementsperproc=0;
154 break;
155 }
156 }
157 }
158 alllastrow[nproc-1]=na-1;
159 for(i=0;i<nproc;i++)columnmap[i]=irow[allfirstrow[i]];
160 }
161 /* MPI_Bcast(columnmap,nproc,MPI_INT,0,MPI_COMM_WORLD); */
162 GA_Brdcst(columnmap, sizeof(int)*nproc, 0);
163 /* MPI_Bcast(allfirstrow,nproc,MPI_INT,0,MPI_COMM_WORLD); */
164 GA_Brdcst(allfirstrow, sizeof(int)*nproc, 0);
165 /* MPI_Bcast(alllastrow,nproc,MPI_INT,0,MPI_COMM_WORLD); */
166 GA_Brdcst(alllastrow, sizeof(int)*nproc, 0);
167 myfirstrow = allfirstrow[me];
168 mylastrow = alllastrow[me];
169
170
171 amat = NGA_Create_irreg(MT_C_DBL, 1, &nz , "A", &nproc,columnmap);
172 if(!amat) GA_Error("create failed: A",nz);
173 if(me==0){
174 lo=0;
175 hi=nz-1;
176 NGA_Put(amat,&lo,&hi,a,&hi);
177 }
178
179 /*now create column matrix */
180 cidx = NGA_Create_irreg(MT_C_INT, 1, &nz , "COLUMN",&nproc,columnmap);
181 if(!cidx) GA_Error("create cidx failed",nz);
182 if(me==0){
183 lo=0;
184 hi=nz-1;
185 NGA_Put(cidx,&lo,&hi,icol,&hi);
186 }
187 GA_Print_distribution(cidx);
188 GA_Sync();
189 /* ************** */
190
191 /*now create row matrix and all vectors*/
192
193 ridx = NGA_Create_irreg(MT_C_INT, 1, &na , "ROW", &nproc,allfirstrow);
194 if(!ridx) GA_Error("create ridx failed",na);
195 lo=0;
196 hi=na-1;
197 if(me==0){
198 NGA_Put(ridx,&lo,&hi,irow,&hi);
199 }
200 GA_Sync();
201 GA_Print_distribution(ridx);
202 MP_BARRIER();
203
204 xvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "X",&nproc,allfirstrow);
205 if(!xvec) GA_Error("create x failed",na);
206 GA_Fill(xvec, &d_one);
207 /*GA_Print(xvec);*/
208
209 bvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "B",&nproc,allfirstrow);
210 if(!bvec) GA_Error("create b failed",na);
211 if(me==0){
212 lo=0;
213 hi=na-1;
214 NGA_Put(bvec,&lo,&hi,x,&hi);
215 }
216 /*GA_Print(xvec);*/
217
218 axvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "Ax", &nproc,allfirstrow);
219 if(!axvec) GA_Error("create Ax failed",na);
220
221 rvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "R", &nproc, allfirstrow);
222 if(!rvec) GA_Error("create r failed",na);
223
224 dvec = GA_Duplicate(rvec,"D"); /* d = r */
225
226 svec = NGA_Create_irreg(MT_C_DBL, 1, &na , "S", &nproc, allfirstrow);
227 if(!svec) GA_Error("create s failed",na);
228
229 dmvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "Dm", &nproc, allfirstrow);
230 if(!dmvec) GA_Error("create dm failed",na);
231
232 #if 0
233 dims[0]=1;dims[1]=na;
234 t_rvec = NGA_Create(MT_C_DBL, 2, dims, "T_R", NULL);
235 if(!t_rvec) GA_Error("create q failed",na);
236 #endif
237
238 qvec = NGA_Create_irreg(MT_C_DBL, 1, &na , "Q", &nproc, allfirstrow);
239 if(!qvec) GA_Error("create q failed",na);
240 /*GA_Print_distribution(qvec);*/
241
242 if(isvectormirrored){
243 ph= GA_Pgroup_get_mirror();
244 m_dvec =
245 NGA_Create_irreg_config(MT_C_DBL,1,&na,"mD",&nproc,allfirstrow,ph);
246 if(!m_dvec) GA_Error("create mirrored dvec failed",na);
247 }
248
249 /* JAD What is GA_Malloc_local? Where is it? */
250 /* ga_vecptr = (double *)GA_Malloc_local(na*sizeof(double)); */
251
252
253 if(me==0)fclose(fd);
254 /*dont forget to free mallocs*/
255 free(allfirstrow);
256 free(alllastrow);
257 free(columnmap);
258 free(a);
259 free(x);
260 free(irow);
261 free(icol);
262 }
263