00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 double arclengthrv(long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp)
00030 {
00031 long i,j,k,n,ni,ini,stop,stopj,modif,li,numr;
00032 double dl, dlmax, dlmin, psi;
00033 double lambda, blambda, dlambda, ddlambda;
00034 double norf, norfp, norv, norfa;
00035 double zero, ierr, alerr;
00036 double dlambdaprev;
00037 double *r;
00038 double *ra;
00039 double *ddr;
00040 double *u;
00041 double *v;
00042 double *f;
00043 double *fa;
00044 double *fp;
00045 double *fc;
00046 double *fi;
00047 const char *mode;
00048 double totl = 0.0;
00049 long back_dl = 0;
00050 long bcid = 0;
00051 long sm_change = 0;
00052 long modif_psi = 0;
00053 flagsw check_rv_fp = nlman->check_lambdar_p;
00054 double bdl;
00055 matrix lsm_a(3,3);
00056 vector lsm_r(3), lsm_l(3);
00057
00058
00059
00060 n = Ndofm;
00061
00062 zero = Mp->zero;
00063
00064
00065 ni = nlman->nial;
00066
00067 ini = nlman->niilal;
00068
00069 ierr = nlman->erral;
00070
00071 dl = nlman->dlal;
00072
00073 alerr = nlman->alerr;
00074
00075 dlmax = nlman->dlmaxal;
00076
00077 dlmin = nlman->dlminal;
00078
00079 psi = nlman->psial;
00080
00081
00082 r = new double[n];
00083 ddr = new double[n];
00084 u = new double[n];
00085 v = new double[n];
00086 f = new double[n];
00087 fi = new double[n];
00088
00089
00090
00091 dlambda=0.0;
00092
00093 if (Mp->hdbcont.restore_stat())
00094 {
00095 if (Mespr==1)
00096 fprintf (stdout,"\n Reading of backup file\n");
00097 solver_restore(ra, fp, li, lambda, dl, NULL, n);
00098
00099 Mm->initmaterialmodels();
00100 print_init(-1, "at");
00101 }
00102 else
00103 {
00104 lambda=0.0;
00105 li=0;
00106 print_init(-1, "wt");
00107 print_step(lcid, li, lambda, f);
00108 }
00109
00110
00111
00112 assemble_attained_load_vector(fa, fc, fp, n, lambda, nlman);
00113
00114 norfp = normv (fp,n);
00115
00116 modif=0;
00117
00118
00119
00120 if (psi < 0) modif_psi = 1;
00121
00122
00123 dlambdaprev=1.0;
00124
00125
00126
00127
00128 for (i=li;i<ni;i++)
00129 {
00130
00131
00132 copyv(ra, r, n);
00133
00134 blambda=lambda;
00135
00136 if (Mespr==1)
00137 {
00138 fprintf (stdout,"\n increment %ld, lambda=%e dl=%e psi=%e",i,lambda,dl,psi);
00139
00140 }
00141
00142
00143 sm_change = assemble_stiffness_matrix (lcid,i,-1,li);
00144
00145
00146 copyv(fp, f, n);
00147
00148
00149 if (sm_change)
00150 Mp->ssle->solve_system(Gtm,Smat,v,f,Out);
00151
00152
00153 norv = displincr(v,n);
00154
00155 if ((i == 0) && (modif_psi) && (norfp > Mp->zero)){
00156
00157 psi0 = psi = sqrt(norv*norv/norfp/norfp);
00158 }
00159
00160
00161 tmp = norv*norv+psi*psi*norfp*norfp;
00162
00163 if (dlambdaprev>0.0){
00164 dlambda = dl/sqrt(tmp);
00165 }else{
00166 dlambda = 0.0-dl/sqrt(tmp);
00167 }
00168
00169
00170 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar_p))
00171 {
00172 check_rv_fp = off;
00173 dlambda = nlman->lambdar_p - lambda;
00174
00175 bdl = dl;
00176
00177 dl = dlambda*sqrt(tmp);
00178 back_dl = 1;
00179 }
00180
00181 if (Mespr == 1)
00182 fprintf(Out, "Step %ld: ddlambda = %e, norv = %e, norfp = %e\n", i, ddlambda, norv, norfp);
00183
00184
00185 copymultv(v, ddr, dlambda, n);
00186
00187
00188 addv(ra, ddr, n);
00189
00190
00191
00192 update_attained_load_vector(fa, fc, fp, n, lambda, dlambda, nlman);
00193
00194
00195 ddlambda=dlambda;
00196
00197
00198 internal_forces (lcid,fi);
00199
00200
00201 subv(fa, fi, f, n);
00202
00203
00204 norf = normv(f,n);
00205 norfa = normv(fa, n);
00206
00207
00208
00209 if (fabs(norfa) > Mp->zero)
00210 norf /= norfa;
00211
00212 if (norf<ierr)
00213 {
00214
00215
00216
00217 if (Mespr==1)
00218 fprintf(stdout,"\n increment %ld norf=%e ierr=%e ddlambda=%e dlambda=%e",i,norf,ierr, ddlambda, dlambda);
00219
00220
00221 totl += dl;
00222
00223
00224
00225 if (back_dl)
00226 dl = bdl;
00227 else
00228 modif++;
00229
00230 if (modif>1)
00231 {
00232
00233 dl*=2.0;
00234 if (dl>dlmax)
00235 dl=dlmax;
00236 modif=0;
00237 if (Mespr==1)
00238 fprintf(stdout,"\n arc-length is be modified (dl increase) dl=%e",dl);
00239 }
00240 lambda+=dlambda;
00241 dlambdaprev=dlambda;
00242
00243
00244 if (modif_psi)
00245 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi);
00246
00247 Mm->updateipval();
00248 compute_req_val(lcid);
00249 print_step(lcid, i, lambda, fa);
00250 print_flush();
00251
00252 if ((nlman->check_lambdar_p == on) && (check_rv_fp == off))
00253 {
00254
00255 break;
00256 }
00257 if ((check_max_tot_al==on) && (fabs(totl - nlman->max_tot_al) < nlman->errtal))
00258 {
00259
00260 break;
00261 }
00262 continue;
00263 }
00264
00265
00266
00267
00268 if (Mespr==1)
00269 fprintf (stdout,"\n inner loop: 1 norf=%e ierr=%e ddlambda=%e dlambda=%e",norf,ierr,ddlambda, dlambda);
00270
00271
00272 stop=0;
00273
00274 stopj=0;
00275
00276 fillm(0.0, lsm_a);
00277 fillv(0.0, lsm_r);
00278
00279 for (j=0;j<ini;j++)
00280 {
00281
00282
00283 Mp->jstep = j;
00284
00285
00286 Mespr = 0;
00287 sm_change = assemble_stiffness_matrix (lcid,i,j,li);
00288 Mespr = 1;
00289
00290
00291 Mp->ssle->solve_system(Gtm,Smat,u,f,Out);
00292 if (sm_changed)
00293 {
00294
00295 copyv(fp, f, n);
00296 Mp->ssle->solve_system(Gtm,Smat,v,f,Out);
00297 }
00298
00299
00300
00301 copyv(ddr, f, n);
00302
00303
00304
00305 addv(ddr, u, n);
00306
00307
00308 determine_dlambda (ddr, u, v, f, n, ddlambda,psi,norfp,dl,dlambda,stopj, nlman);
00309 if (stopj)
00310 {
00311 break;
00312 }
00313
00314
00315 if ((nlman->check_lambdar_p == on) && (check_rv_fp == off) &&
00316 (dlambda+ddlambda+lambda > nlman->lambdar_p))
00317 dlambda = lambdar_p - ddlambda - lambda;
00318
00319
00320
00321 addmultv(ddr, v, dlambda, n);
00322
00323
00324
00325 addv(ra, u, n);
00326 addmultv(ra, v, dlambda, n);
00327
00328
00329 update_attained_load_vector(fa, fc, fp, n, lambda+ddlambda, dlambda, nlman);
00330
00331 ddlambda+=dlambda;
00332
00333
00334
00335 internal_forces (lcid,fi);
00336
00337 subv(fa, fi, f, n);
00338
00339 norf=sizev(f,n);
00340
00341 norfa = sizev(fa, n);
00342 if (fabs(norfa) > zero)
00343 {
00344
00345
00346 norf /= norfa;
00347 }
00348
00349 if (Mp->nlman->dlam = fullmethod){
00350
00351 noral = displincr (ddr,ddr,n);
00352 noral += ddlambda*ddlambda*norfp*norfp*psi*psi;
00353 noral -= dl*dl;
00354 }
00355
00356 if (Mespr==1)
00357 {
00358 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00359
00360 }
00361
00362
00363 if (norf<ierr)
00364 {
00365 dlambdaprev=ddlambda;
00366 lambda+=ddlambda;
00367 compute_req_val (lcid);
00368 print_step(lcid, i, lambda, fa);
00369 print_flush();
00370 Mm->updateipval();
00371 totl += dl;
00372 stop=1;
00373 break;
00374 }
00375
00376
00377 stopj = check_divergence(nlman, lsm_a, lsm_r, lsm_l, j);
00378 if (stopj)
00379 break;
00380
00381 }
00382
00383
00384 modif=0;
00385
00386
00387 if ((check_max_tot_al==on) && (fabs(totl - nlman->max_tot_al) < nlman->errtal))
00388 break;
00389
00390 if (stop==1)
00391 {
00392
00393 if (modif_psi)
00394 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi);
00395
00396 if ((nlman->check_lambdar_p == on) && (check_rv_fp == off) && (fabs(lambda-nlman->lambdar_p) > nlman->errl))
00397 check_rv_fp = on;
00398
00399 if (back_dl)
00400 {
00401 dl = bdl;
00402 back_dl = 0;
00403 }
00404 }
00405
00406 if ((stop==0) || (stopj))
00407
00408 {
00409 if (dl == dlmin)
00410 break;
00411
00412 dl/=2.0;
00413 if (dl<dlmin)
00414 dl=dlmin;
00415
00416 if (Mespr==1)
00417 {
00418 if (dl == dlmin)
00419 {
00420 fprintf (stdout,"\n arc length was decreased dl=dlmin=%e",dl);
00421
00422 }
00423 else
00424 {
00425 fprintf (stdout,"\n arc length was decreased dl=%e",dl);
00426
00427 }
00428 }
00429
00430 copyv (r, ra, n);
00431
00432 lambda=blambda;
00433 }
00434
00435
00436
00437
00438 }
00439
00440 if (Mp->hdbcont.save_stat())
00441 {
00442 if (Mespr==1)
00443 fprintf (stdout,"\n Creating backup file\n");
00444 solver_restore(ra, fp, li, lambda, dl, NULL, n);
00445 }
00446
00447 print_close();
00448
00449 delete [] r;
00450 delete [] fi;
00451 delete [] f;
00452 delete [] v;
00453 delete [] u;
00454 delete [] ddr;
00455
00456 return (lambda);
00457 }
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 long assemble_stiffness_matrix (long lcid,long i,long j,long li)
00481 {
00482 long ret = 0;
00483
00484 switch (Mp->nlman->stmat){
00485 case initial_stiff:{
00486 if (i==li && j<0){
00487
00488 stiffness_matrix (lcid);
00489 ret = 1;
00490 }
00491 break;
00492 }
00493 case tangent_stiff:{
00494 stiffness_matrix (lcid);
00495 ret = 1;
00496 break;
00497 }
00498 case secant_stiff:{
00499 stiffness_matrix (lcid);
00500 ret = 1;
00501 break;
00502 }
00503 case incr_tangent_stiff:{
00504 if (j<0){
00505
00506
00507 stiffness_matrix (lcid);
00508 ret = 1;
00509 }
00510 break;
00511 }
00512 case ijth_tangent_stiff:{
00513 if (j<0){
00514
00515 if ((i%Mp->nlman->ithstiffchange == 0) || (i==li))
00516 {
00517 stiffness_matrix (lcid);
00518 ret = 1;
00519 }
00520 }
00521 else{
00522
00523 if (j%Mp->nlman->jthstiffchange == 0)
00524 {
00525 stiffness_matrix (lcid);
00526 ret = 1;
00527 }
00528 }
00529 break;
00530 }
00531 default:{
00532 print_err("unknown type of stiffness matrix is required",__FILE__,__LINE__,__func__);
00533 }
00534 }
00535 return ret;
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549 void seldofinit ()
00550 {
00551 long i,j,k,l,ndofn;
00552 double x1,x2,y1,y2,z1,z2,length;
00553
00554 switch (Mp->nlman->displnorm){
00555 case alldofs:{ break; }
00556 case seldofs:
00557 case seldofscoord:{
00558 for (i=0;i<Mp->nlman->nsdofal;i++){
00559 j=Mp->nlman->seldofal[i];
00560 Mp->nlman->seldofal[i]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1;
00561 }
00562 break;
00563 }
00564 case selecnodes:{
00565 Mp->nlman->nsdofal=0;
00566 for (i=0;i<Mp->nlman->nsnal;i++){
00567 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00568 }
00569 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
00570 k=0;
00571 for (i=0;i<Mp->nlman->nsnal;i++){
00572 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00573 for (j=0;j<ndofn;j++){
00574 l=Mt->give_dof (Mp->nlman->selnodal[i],j);
00575 if (l>0){
00576 Mp->nlman->seldofal[k]=l-1;
00577 k++;
00578 }
00579 }
00580 }
00581 Mp->nlman->nsdofal=k;
00582 break;
00583 }
00584 case nodesdistincr:{
00585 Mp->nlman->nsdofal=0;
00586 for (i=0;i<Mp->nlman->nsnal;i++){
00587 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00588 }
00589 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
00590 k=0;
00591 for (i=0;i<Mp->nlman->nsnal;i++){
00592 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00593 for (j=0;j<ndofn;j++){
00594 Mp->nlman->seldofal[k]=Mt->give_dof (Mp->nlman->selnodal[i],j);
00595 k++;
00596 }
00597 }
00598 x1=Gtm->gnodes[Mp->nlman->selnodal[0]].x; x2=Gtm->gnodes[Mp->nlman->selnodal[1]].x;
00599 y1=Gtm->gnodes[Mp->nlman->selnodal[0]].y; y2=Gtm->gnodes[Mp->nlman->selnodal[1]].y;
00600 z1=Gtm->gnodes[Mp->nlman->selnodal[0]].z; z2=Gtm->gnodes[Mp->nlman->selnodal[1]].z;
00601 length=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
00602 Mp->nlman->nxal=(x2-x1)/length; Mp->nlman->nyal=(y2-y1)/length; Mp->nlman->nzal=(z2-z1)/length;
00603 break;
00604 }
00605 default:{
00606 fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__);
00607 }
00608 }
00609 }
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633 void assamble_attained_load_vector(double *fa, double *fc, double *fp, long n, double lambda, nonlinman *nlman)
00634 {
00635 double aux = lambda;
00636
00637 if ((nlman->check_lambdar_p == on) && (lambda > nlman->lambdar_p))
00638 aux = nlman->lambdar_p
00639 addmultv(fc, fp, aux, fa, n);
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662 void update_attained_load_vector(double *fa, double *fp, long n, double lambda, double dlambda, nonlinman *nlman)
00663 {
00664 double aux = dlambda;
00665
00666 if ((nlman->check_lambdar_p == on) && (lambda+dlambda > nlman->lambdar_p))
00667 {
00668 if (lambda < nlman->lambdar_p)
00669 aux = nlman->lambdar_p-lambda-dlambda;
00670 else
00671 {
00672
00673
00674 return;
00675 }
00676 }
00677 addmultv(fa, fp, aux, n);
00678 }
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707 long check_divergence(nonlinman *nlman, matrix &lsm_a, vector &lsm_r, vector &lsm_l, long j, double norf)
00708 {
00709 if (nlman->check_div == on)
00710 {
00711 if (j >= nlman->div_min_steps)
00712 {
00713 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, Mp->zero,1) > 0.0)
00714 {
00715 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00716 return(1);
00717 }
00718 }
00719 else
00720 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, Mp->zero,0);
00721 }
00722 return(0);
00723 }
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741 double displincr (double *dr1, double *dr2, long n)
00742 {
00743 long i;
00744 double norm,u1,u2,v1,v2,w1,w2;
00745
00746 switch (Mp->nlman->displnorm){
00747 case alldofs:{
00748 norm = ss (dr1,dr2,n);
00749 break;
00750 }
00751 case seldofs:
00752 case seldofscoord:
00753 case selecnodes:{
00754 norm=0.0;
00755 for (i=0;i<Mp->nlman->nsdofal;i++){
00756 norm+=dr1[Mp->nlman->seldofal[i]]*dr2[Mp->nlman->seldofal[i]];
00757 }
00758 break;
00759 }
00760 case nodesdistincr:{
00761 if (Mp->nlman->probdimal==2){
00762 u1=0.0; u2=0.0; v1=0.0; v2=0.0;
00763 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]-1];
00764 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]-1];
00765 if (Mp->nlman->seldofal[2]>0) u2=dr2[Mp->nlman->seldofal[2]-1];
00766 if (Mp->nlman->seldofal[3]>0) v2=dr2[Mp->nlman->seldofal[3]-1];
00767
00768 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal;
00769 }
00770 if (Mp->nlman->probdimal==3){
00771 u1=0.0; u2=0.0; v1=0.0; v2=0.0; w1=0.0; w2=0.0;
00772 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]-1];
00773 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]-1];
00774 if (Mp->nlman->seldofal[2]>0) w1=dr1[Mp->nlman->seldofal[2]-1];
00775 if (Mp->nlman->seldofal[3]>0) u2=dr2[Mp->nlman->seldofal[3]-1];
00776 if (Mp->nlman->seldofal[4]>0) v2=dr2[Mp->nlman->seldofal[4]-1];
00777 if (Mp->nlman->seldofal[5]>0) w2=dr2[Mp->nlman->seldofal[5]-1];
00778
00779 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal + (w2-w1)*Mp->nlman->nzal;
00780 }
00781 break;
00782 }
00783 default:
00784 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__);
00785 }
00786 return norm;
00787 }
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 void quadeqcoeff (double *ddr,double *v,long n,double ddlambda,double psi,double norfp,double dl,
00811 double &a0,double &a1,double &a2)
00812 {
00813 long i;
00814 double norddr, norv, norddrv;
00815
00816 switch (Mp->nlman->displnorm){
00817 case alldofs:{
00818 case seldofs:
00819 case seldofscoord:
00820 case selecnodes:{
00821 a0 = displincr(ddr, ddr, n);
00822 a1 = 2.0*displincr(ddr, v, n);
00823 a2 = displincr(v, v, n);
00824
00825 a0 += ddlambda*ddlambda*psi*psi*norfp*norfp - dl*dl;
00826 a1 += 2.0*ddlambda*psi*psi*norfp*norfp;
00827 a2 += psi*psi*norfp*norfp;
00828 break;
00829 }
00830 case nodesdistincr:{
00831 norddr = displincr(ddr, ddr, n);
00832 norv = displincr(v, v, n);
00833 a0 = norddr*norddr + ddlambda*ddlambda*norfp*norfp*psi*psi-dl*dl;
00834 a1 = 2.0*(norv*norddr + ddlambda*norfp*norfp*psi*psi);
00835 a2 = norv*norv + norfp*norfp*psi*psi;
00836 break;
00837 }
00838 default:
00839 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__);
00840 }
00841 }
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867 void determine_dlambda (double *ddr, double *u, double *v, double *ddrprev, long n,
00868 double ddlambda, double psi, double norfp, double dl,
00869 double &dlambda, long &stop, nonlinman *nlman)
00870 {
00871 long k,numr;
00872 double l1,l2,aux,a0,a1,a2,ss1,ss2,ss3,ss4,ss5,nom,denom;
00873
00874 switch (Mp->nlman->dlam)
00875 {
00876 case nodetermination:{
00877 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__);
00878 break;
00879 }
00880 case minvalue:
00881 case maxvalue:
00882 case minangle:{
00883
00884 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
00885
00886 numr = solv_polynom_2(a2, a1, a0, l1, l2);
00887 switch (numr){
00888 case -1:{
00889 print_err("infinite number of solution of constrained condition\n"
00890 "(all coefficients of quadratic equation are equal to zero", __FILE__, __LINE__, __func__);
00891 abort();
00892 break;
00893 }
00894 case 0:{
00895 print_err("nonsolvable constrained condition in function arclength", __FILE__, __LINE__, __func__);
00896 stop=0;
00897 break;
00898 }
00899 case 1:{
00900 dlambda = l1;
00901 break;
00902 }
00903 default:{}
00904 }
00905 break;
00906 }
00907 case linearizedmeth:{
00908 break;
00909 }
00910 case fullmethod:
00911 break;
00912 default:{
00913 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__);
00914 }
00915 }
00916
00917
00918 switch (Mp->nlman->dlam)
00919 {
00920 case nodetermination:{
00921 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__);
00922 break;
00923 }
00924 case minvalue:{
00925 if (fabs(l1)>fabs(l2))
00926 dlambda=l2;
00927 else
00928 dlambda=l1;
00929 break;
00930 }
00931 case maxvalue:{
00932 if (fabs(l1)>fabs(l2))
00933 dlambda=l1;
00934 else
00935 dlambda=l2;
00936 break;
00937 }
00938 case minangle:{
00939 ss1=0.0; ss2=0.0;
00940 ss3=0.0; ss4=0.0; ss5=0.0;
00941 for (k=0;k<n;k++){
00942 ss1+=(ddr[k]+l1*v[k])*ddrprev[k];
00943 ss2+=(ddr[k]+l2*v[k])*ddrprev[k];
00944 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00945 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00946 ss5+=ddrprev[k]*ddrprev[k];
00947 }
00948 if (ss1/sqrt(ss3)/sqrt(ss5)>ss2/sqrt(ss4)/sqrt(ss5))
00949 dlambda=l1;
00950 else
00951 dlambda=l2;
00952 break;
00953 }
00954 case linearizedmeth:{
00955
00956 lij = displincr (ddrprev,ddrprev,n);
00957 lij += psi*psi*norfp*norfp*ddlambda;
00958 lij -= dl*dl;
00959
00960 nom = -0.5*lij - displincr(ddrprev,u,n);
00961 denom = displincr(ddrprev,v,n);
00962 denom += psi*psi*norfp*norfp*ddlambda;
00963
00964 dlambda=nom/denom;
00965 break;
00966 }
00967 case fullmethod:{
00968
00969 lij = displincr (ddrprev,ddrprev,n);
00970 lij += psi*psi*norfp*norfp*ddlambda;
00971 lij -= dl*dl;
00972
00973
00974
00975
00976
00977
00978 nom = 0.0 - displincr (ddrprev,u,n);
00979 nom -= lij;
00980
00981 denom = 2.0*displincr (ddrperv,v,n);
00982 denom += 2.0*psi*psi*ddlambda*norfp*norfp;
00983
00984 dlambda = nom/denom;
00985 break;
00986 }
00987 default:{
00988 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__);
00989 }
00990 }
00991
00992 }
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016 void adapt_psi_coeff(long i, double ddlambda, double &ddlambda0, double &psi0, double &psi)
01017 {
01018
01019
01020 if (i==0)
01021 ddlambda0 = ddlambda;
01022
01023
01024
01025 if ((sgn(ddlambda0) != sgn(ddlambda)) && (sgn(ddlambda) != 0.0))
01026 {
01027 psi0= psi;
01028 ddlambda0 = ddlambda;
01029 }
01030
01031
01032
01033 if (ddlambda0 != 0.0)
01034 psi = psi0*fabs(ddlambda/ddlambda0);
01035
01036 return;
01037 }