00001 #include <errno.h>
00002
00003 #include "slesolv.h"
00004 #include "mathem.h"
00005
00006 #include "iotools.h"
00007
00008 #include "precond.h"
00009 #include "gmatrix.h"
00010 #include "gtopology.h"
00011 #include "seqselnodes.h"
00012 #include "saddpoint.h"
00013 #include "aggregator.h"
00014
00015 slesolv::slesolv ()
00016 {
00017
00018 tlinsol=gauss_elim;
00019 tsol=gauss_elim;
00020
00021 zero=1.0e-15;
00022
00023
00024 ni=0;
00025
00026 err=0.0;
00027
00028 ani=0;
00029
00030 aerr=0.0;
00031
00032
00033
00034 iv = 0;
00035
00036
00037 nbdof=0;
00038
00039
00040
00041 bsize=0;
00042
00043 }
00044
00045 slesolv::~slesolv ()
00046 {
00047
00048 }
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 void slesolv::read (gtopology *top,XFILE *in,long mespr)
00060 {
00061
00062 xfscanf (in,"%k%m","typelinsol",&linsolvertype_kwdset,(int*)&tlinsol);
00063
00064 switch (tlinsol){
00065 case gauss_elim:{
00066 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the gaussian elimination method");
00067 break;
00068 }
00069 case ldl:{
00070 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the LDL method");
00071 break;
00072 }
00073 case lu:{
00074 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the LU method");
00075 break;
00076 }
00077 case ll:{
00078 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the LL method (Cholesky)");
00079 break;
00080 }
00081 case ill:{
00082 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the incomplete LL method (Cholesky)");
00083 break;
00084 }
00085 case conden:{
00086 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by condensation method");
00087 xfscanf (in,"%k%m","typesecondsol",&linsolvertype_kwdset,(int*)&tsol);
00088 xfscanf (in,"%k%ld","nbdof",&nbdof);
00089 break;
00090 }
00091 case cg:{
00092 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the CG method");
00093
00094
00095 xfscanf (in,"%k%ld","number_of_iterations",&ni);
00096 xfscanf (in,"%k%lf","error_of_computation",&err);
00097 xfscanf (in,"%k%lf","initial_values",&iv);
00098 break;
00099 }
00100 case bicg:{
00101 if (mespr==1) fprintf (stdout,"\n system of linear equationes will be solved by the bi-conjugate gradient method");
00102
00103
00104 xfscanf (in,"%k%ld","number_of_iterations",&ni);
00105 xfscanf (in,"%k%lf","error_of_computation",&err);
00106 xfscanf (in,"%k%lf","initial_values",&iv);
00107 break;
00108 }
00109
00110 case lapack_sol:{
00111 break;
00112 }
00113
00114 case spdirldl:{
00115 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by sparse direct method (LDL)");
00116 break;
00117 }
00118 case spdirlu:{
00119 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by sparse direct method (LU)");
00120 break;
00121
00122 }
00123 case spdirll:{
00124 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by sparse direct method (LL)");
00125 break;
00126 }
00127
00128 case pardisolver:{
00129
00130 print_err("PARDISO solver is not supported at this time",__FILE__,__LINE__,__func__);
00131 break;
00132 }
00133
00134 case sschur:{
00135 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the Schur complement method");
00136
00137 schur.read (top,mespr,in);
00138
00139 top->rst=1;
00140 top->cngen=0;
00141
00142 break;
00143 }
00144 case sfeti:{
00145 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the FETI method");
00146
00147 feti.read (top,in);
00148
00149 break;
00150 }
00151 case saddle_point:{
00152 if (mespr==1) fprintf (stdout,"\n system of linear equations will be solved by the saddle point method");
00153
00154
00155
00156
00157
00158 break;
00159 }
00160 default:{
00161 print_err("unknown type of solver of linear equations is required",__FILE__,__LINE__,__func__);
00162 }
00163 }
00164
00165 fprintf (stdout,"\n computer zero in slesolv is defined as %le",zero);
00166
00167 switch (tlinsol){
00168 case gauss_elim:
00169 case ldl:
00170 case lu:
00171 case ll:
00172 case ill:
00173 case conden:
00174 case lapack_sol:
00175 case spdirldl:
00176 case spdirlu:
00177 case spdirll:
00178 case pardisolver:{
00179
00180 break;
00181 }
00182 case cg:
00183 case bicg:{
00184 prec.read (top,in,mespr);
00185 break;
00186 }
00187 case saddle_point:
00188 case sschur:
00189 case sfeti:{
00190 break;
00191 }
00192 default:{
00193 print_err("unknown type of solver of linear equations is required",__FILE__,__LINE__,__func__);
00194 }
00195 }
00196
00197
00198 if (prec.pt==boss){
00199 if (prec.agg->impl==2){
00200
00201 top->rst=1;
00202 }
00203 }
00204
00205
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215 void slesolv::print (FILE *out)
00216 {
00217 fprintf (out,"%d\n",(int)tlinsol);
00218
00219 switch (tlinsol){
00220 case gauss_elim:
00221 case ldl:
00222 case lu:
00223 case ll:
00224 case ill:
00225 case lapack_sol:
00226 case spdirldl:
00227 case spdirlu:
00228 case spdirll:
00229 case pardisolver:{
00230 break;
00231 }
00232 case conden:{
00233 fprintf (out,"%ld\n",nbdof);
00234 break;
00235 }
00236 case sfeti:{
00237 feti.print (out);
00238 break;
00239 }
00240
00241 case cg:{
00242 fprintf (out,"%ld %e %ld\n", ni,err,iv);
00243 break;
00244 }
00245 case bicg:{
00246 fprintf (out,"%ld %e %ld\n",ni,err,iv);
00247 break;
00248 }
00249 case saddle_point:{
00250 break;
00251 }
00252 default:{
00253 print_err("unknown type of solver of linear equations is required",__FILE__,__LINE__,__func__);
00254 }
00255 }
00256
00257
00258 switch (tlinsol){
00259 case gauss_elim:
00260 case ldl:
00261 case lu:
00262 case ll:
00263 case ill:
00264 case lapack_sol:
00265 case spdirldl:
00266 case spdirlu:
00267 case spdirll:
00268 case pardisolver:{
00269
00270 break;
00271 }
00272 case cg:
00273 case bicg:{
00274 prec.print (out);
00275 break;
00276 }
00277 case sfeti:
00278 case saddle_point:{
00279 break;
00280 }
00281 default:{
00282 print_err("unknown type of solver of linear equations is required",__FILE__,__LINE__,__func__);
00283 }
00284 }
00285
00286 fprintf (out,"\n");
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296 void slesolv::initiate (slesolv *ssle)
00297 {
00298 tlinsol = ssle->tlinsol;
00299 tsol = ssle->tsol;
00300
00301 bsize = ssle->bsize;
00302 zero = ssle->zero;
00303
00304 ni = ssle->ni;
00305 err = ssle->err;
00306 ani = ssle->ani;
00307 aerr = ssle->aerr;
00308 iv = ssle->iv;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 void slesolv::prepare_data (long *ndof,gtopology *top,FILE *out)
00321 {
00322 long i,j;
00323
00324 if (tlinsol == sschur || tlinsol == sfeti){
00325 if (top->stop->md==glob_glued){
00326 top->auto_subdom_detection (out);
00327 }
00328 }
00329
00330
00331
00332
00333 if (tlinsol == sschur){
00334 top->stop->coarse_local_map (top,out);
00335
00336 long **schurnodes;
00337 schurnodes = new long* [top->stop->ns];
00338 for (i=0;i<top->stop->ns;i++){
00339 schurnodes[i] = new long [top->stop->nnsd[i]];
00340 }
00341
00342
00343 for (i=0;i<top->stop->ns;i++){
00344 for (j=0;j<top->stop->nnsd[i];j++){
00345 if (top->stop->nodmultip[i][j]>1)
00346 schurnodes[i][j]=1;
00347 else
00348 schurnodes[i][j]=-1;
00349 }
00350 }
00351
00352 seqselnodes *selnodschur;
00353 selnodschur = new seqselnodes (top->stop->ns,top->stop->nnsd,schurnodes,
00354 top->stop->nodmultip,top->stop->ggnbn,top->stop->icnbnmas,
00355 top->stop->tnbn,top->stop->icmultip,top->stop->lnbncn,top->stop->ggnbncn,top->stop->sid,
00356 top->stop->md,out,1);
00357
00358
00359
00360
00361 selnodschur->prepare_schur (top,out);
00362 ndof[0]=top->stop->schur_ordering (top,out);
00363
00364 schur.initiate (selnodschur);
00365 schur.assemble_subdom_unknowns (top,out);
00366
00367 delete selnodschur;
00368
00369 for (i=0;i<top->stop->ns;i++){
00370 delete [] schurnodes[i];
00371 }
00372 delete [] schurnodes;
00373
00374 }
00375
00376
00377
00378
00379 if (tlinsol == sfeti){
00380 top->stop->coarse_local_map (top,out);
00381
00382 long **fetinodes;
00383 fetinodes = new long* [top->stop->ns];
00384 for (i=0;i<top->stop->ns;i++){
00385 fetinodes[i] = new long [top->stop->nnsd[i]];
00386 }
00387
00388 for (i=0;i<top->stop->ns;i++){
00389 for (j=0;j<top->stop->nnsd[i];j++){
00390 if (top->stop->nodmultip[i][j]>1)
00391 fetinodes[i][j]=1;
00392 else
00393 fetinodes[i][j]=-1;
00394 }
00395 }
00396
00397 seqselnodes *selnodfeti;
00398
00399 selnodfeti = new seqselnodes (top->stop->ns,top->stop->nnsd,fetinodes,
00400 top->stop->nodmultip,top->stop->ggnbn,top->stop->icnbnmas,
00401 top->stop->tnbn,top->stop->icmultip,top->stop->lnbncn,top->stop->ggnbncn,top->stop->sid,
00402 top->stop->md,out,1);
00403
00404 ndof[0]=selnodfeti->prepare_feti (feti.fetiimpl,top,out);
00405
00406 feti.initiate (selnodfeti,top,out);
00407
00408 feti.assemble_subdom_unknowns (top,out);
00409
00410
00411
00412 for (i=0;i<top->stop->ns;i++){
00413 delete [] fetinodes[i];
00414 }
00415 delete [] fetinodes;
00416 }
00417
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430 void slesolv::solve_system (gtopology *gt,gmatrix *gm,double *lhs,double *rhs,FILE *out)
00431 {
00432 time_t bt,et;
00433 long hod,min;
00434 double sec = clock();
00435
00436
00437
00438 prec.initiation (gt,gm,out);
00439
00440 switch (tlinsol){
00441 case gauss_elim:
00442 case ldl:
00443 case lu:
00444 case ll:
00445 case ill:
00446 case cg:
00447 case bicg:
00448 case lapack_sol:
00449 case spdirldl:
00450 case spdirlu:
00451 case spdirll:
00452 case pardisolver:{
00453 bt = time (NULL);
00454
00455 gm->solve_system (gt,prec,lhs,rhs,out);
00456
00457 et = time (NULL);
00458 sec = (clock() - sec) / (double)CLOCKS_PER_SEC;
00459 hod = (long)sec/3600; sec -= hod*3600;
00460 min = (long)sec/60; sec -= min*60;
00461
00462
00463
00464
00465 break;
00466 }
00467 case conden:{
00468 gt->gnodes[5].ai=0;
00469 gt->nidof=gm->n-nbdof;
00470 gm->nbdof=nbdof;
00471 gm->solve_system (gt,prec,lhs,rhs,out);
00472 break;
00473 }
00474 case sschur:{
00475
00476 schur.solve_system (gt,gm,lhs,rhs,out);
00477
00478 break;
00479 }
00480 case sfeti:{
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 feti.solve_system (gt,gm,lhs,rhs,out);
00531
00532 break;
00533 }
00534
00535 case saddle_point:{
00536
00537
00538 break;
00539 }
00540
00541 default:{
00542 print_err("unknown type of solver of linear equations is required",__FILE__,__LINE__,__func__);
00543 }
00544 }
00545 check_math_err();
00546 }
00547