00001 #include "parcongrad.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005
00006
00007
00008
00009 parcongrad::parcongrad(int np,int mr,long nd, long mes)
00010 {
00011 mespr = mes;
00012 nproc=np;
00013 myrank=mr;
00014 ndom=nd;
00015 ndof=0;
00016 nbdof=0;
00017 nbn=0;
00018 maxnbdof=0;
00019 }
00020
00021 parcongrad::~parcongrad()
00022 {
00023
00024 }
00025
00026
00027 void parcongrad::initiate(partop *ptop,gtopology *top,FILE *out)
00028 {
00029 ndof = top->codenum_generation (out);
00030
00031 maxndof = ptop->maxndof;
00032
00033 tndof=ptop->tndof;
00034
00035 bcn = ptop->gcnd;
00036
00037 ndofdom = ptop->nud;
00038
00039
00040 }
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 void parcongrad::buff_coarse (double *cv,double *buff,long nd)
00053 {
00054 long i,j;
00055
00056
00057 for (i=0;i<ndofdom[nd];i++){
00058 j=bcn[nd][i]-1;
00059 cv[j]+=buff[i];
00060 }
00061 }
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 void parcongrad::coarse_buff (double *cv,double *buff,long nd)
00073 {
00074 long i,j;
00075
00076
00077 for (i=0;i<ndofdom[nd];i++){
00078 j=bcn[nd][i]-1;
00079 buff[i]=cv[j];
00080 }
00081 }
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 void parcongrad::solve_system (partop *ptop,gtopology *top,gmatrix *gm,
00102 long *domproc,double *lhs,double *rhs,FILE *out,long iv)
00103 {
00104 long i,j,k,buffsize;
00105 double nom,denom,alpha,beta;
00106 double *d,*p,*r,*dd,*dp,*dr,*dx,*buff,*master_rhs;
00107 MPI_Status stat;
00108
00109 buffsize=maxndof+1;
00110
00111
00112 p = new double [ndof];
00113
00114 r = new double [ndof];
00115
00116 d = new double [ndof];
00117
00118 buff = new double [buffsize];
00119
00120
00121
00122 if (myrank==0){
00123 for (i=1;i<nproc;i++){
00124 MPI_Send (&ni,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00125 }
00126 }
00127 else{
00128 MPI_Recv (&ni,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00129 }
00130 MPI_Barrier (MPI_COMM_WORLD);
00131
00132
00133
00134 if(myrank == 0){
00135 master_rhs = new double[tndof];
00136 nullv (master_rhs,tndof);
00137
00138 k=domproc[0];
00139 buff_coarse (master_rhs,buff,k);
00140
00141 for (j=1;j<nproc;j++){
00142 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00143 k=domproc[stat.MPI_TAG];
00144 buff_coarse (master_rhs,buff,k);
00145 }
00146
00147 if(mespr == 1) fprintf (out,"\n\n\n\n Check of master RHS\n\n");
00148 double norm=0.0;
00149 for(i = 0; i < tndof; i++){
00150 norm+=master_rhs[i];
00151 }
00152 if(mespr == 1){
00153 for(i = 0; i < tndof; i++){
00154 fprintf (out,"%ld %le\n",i,master_rhs[i]);
00155 }
00156 }
00157 if(mespr == 1) fprintf (stdout,"\n\n RHS size is %le",norm);
00158
00159 delete []master_rhs;
00160 }
00161 else{
00162 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00163 }
00164 MPI_Barrier (MPI_COMM_WORLD);
00165
00166
00167
00168
00169
00170
00171
00172 if (iv==0){
00173 for (i=0;i<ndof;i++){
00174 lhs[i]=0.0;
00175 }
00176 }
00177
00178
00179 gm->gmxv (lhs,p);
00180
00181
00182 subv (rhs,p,r,ndof);
00183
00184 nullv (buff,buffsize);
00185 copyv (r,buff,ndof);
00186
00187 if (myrank==0){
00188
00189
00190
00191 dr = new double [tndof];
00192 nullv (dr,tndof);
00193
00194 dd = new double [tndof];
00195 nullv (dd,tndof);
00196
00197 dx = new double [tndof];
00198 nullv (dx,tndof);
00199
00200 dp = new double [tndof];
00201 nullv (dp,tndof);
00202
00203
00204 k=domproc[0];
00205 buff_coarse (dr,buff,k);
00206
00207
00208 for (i=1;i<nproc;i++){
00209 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00210 k=domproc[stat.MPI_TAG];
00211 buff_coarse (dr,buff,k);
00212 }
00213
00214 nom = ss (dr,dr,tndof);
00215
00216 copyv (dr,dd,tndof);
00217
00218
00219 for (i=1;i<nproc;i++){
00220 k=domproc[i];
00221 nullv (buff,buffsize);
00222 coarse_buff (dd,buff,k);
00223 MPI_Send (buff,buffsize,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
00224 }
00225
00226 k=domproc[0];
00227 nullv (buff,buffsize);
00228 coarse_buff (dd,buff,k);
00229 }
00230 else{
00231 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00232 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00233 }
00234 MPI_Barrier (MPI_COMM_WORLD);
00235
00236
00237 copyv (buff,d,ndof);
00238
00239
00240 for (i=0;i<ni;i++){
00241
00242
00243 gm->gmxv (d,p);
00244
00245 copyv (p,buff,ndof);
00246
00247
00248 if (myrank==0){
00249 nullv (dp,tndof);
00250
00251
00252 k=domproc[0];
00253 buff_coarse (dp,buff,k);
00254
00255
00256 for (j=1;j<nproc;j++){
00257 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00258 k=domproc[stat.MPI_TAG];
00259 buff_coarse (dp,buff,k);
00260 }
00261
00262 denom = ss (dd,dp,tndof);
00263
00264 alpha = nom/denom;
00265
00266
00267
00268 for (j=0;j<tndof;j++){
00269 dx[j]+=alpha*dd[j];
00270 dr[j]-=alpha*dp[j];
00271 }
00272
00273 denom=nom;
00274
00275 nom = ss (dr,dr,tndof);
00276
00277 if(mespr == 1) fprintf (stdout,"\n iteration number %ld nom %le",i,nom);
00278
00279 if (nom<err)
00280 buff[buffsize-1]=100.0;
00281
00282 beta = nom/denom;
00283
00284
00285 for (j=0;j<tndof;j++){
00286 dd[j]=beta*dd[j]+dr[j];
00287 }
00288
00289
00290 for (j=1;j<nproc;j++){
00291 k=domproc[j];
00292 nullv (buff,buffsize-1);
00293 coarse_buff (dd,buff,k);
00294 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00295 }
00296
00297 k=domproc[0];
00298 nullv (buff,buffsize-1);
00299 coarse_buff (dd,buff,k);
00300
00301 }
00302 else{
00303 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00304 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00305 }
00306 MPI_Barrier (MPI_COMM_WORLD);
00307
00308
00309 copyv (buff,d,ndof);
00310
00311 if (buff[buffsize-1]>1.0){
00312
00313 if (myrank==0){
00314 for (j=1;j<nproc;j++){
00315 k=domproc[j];
00316 nullv (buff,buffsize-1);
00317 coarse_buff (dx,buff,k);
00318 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00319 }
00320 k=domproc[0];
00321 nullv (buff,buffsize-1);
00322 coarse_buff (dx,buff,k);
00323 }
00324 else{
00325 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00326 }
00327
00328 copyv (buff,lhs,ndof);
00329
00330 break;
00331 }
00332 }
00333 MPI_Barrier (MPI_COMM_WORLD);
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 }