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 "armci.h"
28 #include "message.h"
29
30 extern int na;
31 extern int nz;
32 extern double *dmvec,*svec,*bvec,*dvec,*amat,*xvec,*axvec,*rvec,*qvec;
33 extern int *ridx,*cidx;
34 extern int me, nproc;
35 extern int myfirstrow,mylastrow;
36 static int *columnmap,*allfirstrow,*alllastrow;
37 static FILE *fd;
38
generate_random_file(int naa,int nnz)39 void generate_random_file(int naa,int nnz){
40 fd = fopen("randominput.dat", "w");
41 }
42
read_and_create(int argc,char ** argv)43 void read_and_create(int argc, char **argv)
44 {
45 int ri,i;
46 int tmp1,idealelementsperproc;
47 void **amatptrs,**xvecptrs;
48
49 na = atoi(argv[1]);
50 nz = atoi(argv[2]);
51
52 if(strncmp("random",argv[3],6)){
53 if(me==0){
54 fd = fopen(argv[3], "r");
55 if(fd==NULL)ARMCI_Error("unable to open given file",0);
56 }
57 }
58 else{
59 if(na==0 || nz==0){
60 printf("\nERROR:exiting-no input file given and na or nz is 0");
61 fflush(stdout);
62 ARMCI_Finalize();
63 armci_msg_finalize();
64 return;
65 }
66 if(me==0){
67 generate_random_file(na,nz);
68 fd = fopen("randominput.dat", "r");
69 }
70 }
71 if(me==0){
72 if(na==0)
73 fread(&na, sizeof(na), 1, fd);
74 if(nz==0)
75 fread(&nz, sizeof(nz), 1, fd);
76 printf("\nReading CG input\n");
77 printf("Number of rows: %d\n", na);
78 printf("Number of non-zeros: %d\n", nz);
79 }
80
81 armci_msg_bcast(&nz,sizeof(int),0);
82 armci_msg_bcast(&na,sizeof(int),0);
83 armci_msg_barrier();
84
85 amatptrs = (void **)malloc(sizeof(void *)*nproc);
86 xvecptrs = (void **)malloc(sizeof(void *)*nproc);
87 if(xvecptrs==NULL || amatptrs==NULL)
88 ARMCI_Error("xvecptrs amatptrs malloc failed",sizeof(void *)*nproc);
89
90 if(ARMCI_Malloc(amatptrs,((me==0)?(sizeof(double)*nz):0)))
91 ARMCI_Error("amat malloc failed",sizeof(double)*nz);
92 amat = (double *)amatptrs[0];
93
94 if(ARMCI_Malloc(amatptrs,((me==0)?(sizeof(int)*(nz+1)):0)))
95 ARMCI_Error("icol malloc failed",sizeof(int)*(nz+1));
96 cidx = (int *)amatptrs[0];
97
98 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(int)*(na+1)):0)); /*+1 for end of last row*/
99 ridx = (int *)xvecptrs[0];
100
101 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*(na+1)):0));
102 xvec = (double *)xvecptrs[0];
103
104 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*(na+1)):0));
105 bvec = (double *)xvecptrs[0];
106
107 if(me==0){
108
109 for (i = 0; i < na + 1; i++)
110 xvec[i] = 0.0;
111
112 fread(amat, sizeof(double), nz, fd);
113 fread(ridx, sizeof(int), (na+1), fd);
114 ridx[na]=nz;
115 fread(cidx, sizeof(int), (nz+1), fd);
116 fread(bvec, sizeof(double), (na+1), fd);
117
118 /* the c adjustment */
119 for (i = 0; i < na; i++)
120 ridx[i] -= 1;
121
122 for (i = 0; i < nz; i++)
123 cidx[i] -= 1;
124 }
125
126 armci_msg_barrier();
127 /*acg_matvecmul(amat,xvec,bvec,ridx,cidx);*/
128 if(0){
129 for(i=0;i<nz+1;i++)
130 printf("\n%d:amat[%d]=%f icol[%d]=%d",me,i,amat[i],i,cidx[i]);
131 for(i=0;i<na+1;i++)
132 printf("\n%d:irow[%d]=%d bvec[%d]=%f",me,i,ridx[i],i,bvec[i]);
133 }
134 allfirstrow = (int *)malloc(sizeof(int)*nproc);
135 alllastrow = (int *)malloc(sizeof(int)*nproc);
136 columnmap = (int *)malloc(sizeof(int)*nproc);
137 if(!allfirstrow || !alllastrow || !columnmap)
138 ARMCI_Error("malloc failed allfirstrow ",0);
139 armci_msg_barrier();
140 /*
141 * next decide who works on which rows, this will decide the
142 * distribution of a,d,r,q,x,and ax
143 */
144 /*create the mapping for all vectors, row matrix and column matrix*/
145 if(me==0){
146 idealelementsperproc = nz/nproc;
147 tmp1=0;
148 for(i=0;i<nproc;i++){
149 int elementsperproc=0;
150 allfirstrow[i]=tmp1;
151 for(ri=tmp1;ri<na;ri++,tmp1++){
152 elementsperproc+=(ridx[ri+1]-ridx[ri]);
153 if(elementsperproc>=idealelementsperproc){
154 if((elementsperproc-idealelementsperproc) >
155 idealelementsperproc-(elementsperproc-(ridx[ri+1]-ridx[ri]))){
156 alllastrow[i] = ri-1;
157 if((ri-1)<0)ARMCI_Error("run on a smaller processor count",0);
158 /*tmp1--;*/
159 }
160 else{
161 alllastrow[i] = ri;
162 if(ri<0)ARMCI_Error("run on a smaller processor count",0);
163 tmp1++;
164 }
165 elementsperproc=0;
166 break;
167 }
168 }
169 }
170 alllastrow[nproc-1]=na-1;
171 for(i=0;i<nproc;i++)columnmap[i]=ridx[allfirstrow[i]];
172 }
173 armci_msg_bcast(columnmap,nproc*sizeof(int),0);
174 armci_msg_bcast(allfirstrow,nproc*sizeof(int),0);
175 armci_msg_bcast(alllastrow,nproc*sizeof(int),0);
176 myfirstrow = allfirstrow[me];
177 mylastrow = alllastrow[me];
178 if(me==0)for(i=0;i<nproc;i++){
179 printf("\nDISTRIBUTION:first row of process\t%d is %d last row of process\t%d is %d",i,allfirstrow[i],i,alllastrow[i]);
180 }
181 /*
182 for(i=myfirstrow;i<mylastrow;i++){
183 xvec[i]=0.0;
184 }
185 */
186 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
187 rvec = (double *)xvecptrs[0];
188
189 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
190 dvec = (double *)xvecptrs[0];
191
192 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
193 svec = (double *)xvecptrs[0];
194
195 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
196 dmvec = (double *)xvecptrs[0];
197
198 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
199 qvec = (double *)xvecptrs[0];
200
201 ARMCI_Malloc(xvecptrs,((me==0)?(sizeof(double)*na):0));
202 axvec = (double *)xvecptrs[0];
203
204 if(me==0)fclose(fd);
205 /*dont forget to free mallocs*/
206 free(allfirstrow);
207 free(alllastrow);
208 free(columnmap);
209 }
210