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