1 //#**************************************************************
2 //#
3 //# filename: femath.cpp
4 //#
5 //# author: Gerstmayr Johannes
6 //#
7 //# generated: 15.05.97
8 //# description: Classes for linear and nonlinear algebra which is
9 //# thought to be used in finite element or similar applications
10 //# There are 2D and 3D and arbitrary size Vectors (Vector2D, Vector3D, Vector),
11 //# arbitrary size matrices (Matrix), and a nonlinear system (NumNLSys)
12 //# and a nonlinear solver (NumNLSolver)
13 //# remarks: Indizes run from 1 to n except in Vector3D/2D
14 //#
15 //# Copyright (c) 2003-2013 Johannes Gerstmayr, Linz Center of Mechatronics GmbH, Austrian
16 //# Center of Competence in Mechatronics GmbH, Institute of Technical Mechanics at the
17 //# Johannes Kepler Universitaet Linz, Austria. All rights reserved.
18 //#
19 //# This file is part of HotInt.
20 //# HotInt is free software: you can redistribute it and/or modify it under the terms of
21 //# the HOTINT license. See folder 'licenses' for more details.
22 //#
23 //# bug reports are welcome!!!
24 //# WWW: www.hotint.org
25 //# email: bug_reports@hotint.org or support@hotint.org
26 //#***************************************************************************************
27
28 #include "ioincludes.h"
29
30 #include <assert.h>
31 #include <memory.h>
32
33 #include <math.h>
34
35 #include "tarray.h"
36 #include "mystring.h"
37 #include "femath.h"
38 #include "femathHelperFunctions.h"
39
40
41
42
43 //for output possibility within femath
44 #include "..\workingmodule\WorkingModuleBaseClass.h"
45 extern UserOutputInterface * global_uo;
46 //#include "..\SuperLU\slu_ddefs.h"
47
48 #ifdef gencnt
49 extern int * genvec;
50 extern int * genmat;
51 #endif
52
53 //extern double global_time;
54
55 //ofstream testout("..\\..\\output\\test.txt");
56
57 TArray<double> TMtspent(TMlen);
58 TArray<double> TMtstart(TMlen);
59 TArray<mystr> TMmsg(TMlen);
60
TMResetTimer()61 void TMResetTimer()
62 {
63 TMmsg(1) = mystr("total time");
64 TMmsg(2) = mystr("DrawSystem()");
65 TMmsg(3) = mystr("Full Jacobian");
66 TMmsg(4) = mystr("Partial Jacobian Additioning");
67 TMmsg(5) = mystr("mbs->EvalF2");
68 TMmsg(6) = mystr("mbs->EvalM");
69 TMmsg(7) = mystr("mbs->EvalF");
70 TMmsg(8) = mystr("mbs->EvalG");
71 TMmsg(9) = mystr("Dependencies");
72 TMmsg(10) = mystr("Newton Solve");
73 TMmsg(11) = mystr("Solve");
74 TMmsg(12) = mystr("Mult(Matrix,Vector,Vector&) general");
75 TMmsg(13) = mystr("NLF");
76 TMmsg(14) = mystr("Factorize (sparse+dense jacobi)");
77 TMmsg(15) = mystr("Apply (sparse+dense jacobi)");
78 TMmsg(16) = mystr("dCdq+AddElementLoad (Element::EvalF2)");
79 TMmsg(17) = mystr("Quadratic VV");
80 TMmsg(18) = mystr("Rigid3D-RotMat");
81 TMmsg(19) = mystr("Partial Jacobian");
82 TMmsg(20) = mystr("PostNewtonStep"); //$ LA 2011-2-1: Renameing of old NonlinStep
83 TMmsg(21) = mystr("TimeStep");
84 TMmsg(22) = mystr("EvalF2 CMS");
85 TMmsg(23) = mystr("EvalM CMS");
86 TMmsg(24) = mystr("Build contact search trees");
87 TMmsg(25) = mystr("Contact Search");
88 TMmsg(26) = mystr("Compute gap");
89 TMmsg(27) = mystr("Nearest Segment");
90 TMmsg(28) = mystr("test3");
91 TMmsg(29) = mystr("test4");
92
93 for(int i=1; i<= TMlen; i++)
94 {
95 TMtstart(i) = 0;
96 TMtspent(i) = 0;
97 //(*global_uo) << i << " " << TMmsg(i) << "\n";
98 }
99 }
100
101 //ofstream timings("..\\..\\output\\timings.txt"); //�
102 //ofstream timings2("..\\..\\output\\timings2.txt"); //�
103
104
TMPrintTimer()105 void TMPrintTimer()
106 {
107 #ifdef timeron
108 global_uo->SaveLocalMessageLevel();
109 global_uo->SetLocalMessageLevel(UO_LVL_all);
110
111 double div = TMGetTimer(1);
112 if (TMGetTimer(1) == 0) {div = 1;}
113 for(int i=1; i<= TMlen; i++)
114 {
115 double tt = TMGetTimer(i);
116 if (i>1) {tt=tt/div;}
117 char tstr[30];
118 if (i!=1)
119 {
120 if (tt*100 < 10)
121 sprintf_s(tstr," %06.4f",tt*100);
122 else
123 sprintf_s(tstr,"%06.4f",tt*100);
124 (*global_uo) << mystr("[")+mystr(i)+mystr("]")+mystr(" :\t")+mystr(tstr)+mystr("% : ")+mystr(TMmsg(i))+mystr("\n");
125 }
126 else
127 {
128 sprintf_s(tstr,"%4.2f s",tt);
129 (*global_uo) << mystr("[")+mystr(i)+mystr("]")+mystr(" :\t")+mystr(tstr)+mystr(" : ")+mystr(TMmsg(i))+mystr("\n");
130 }
131 }
132 double acc = TMGetTimer(1);
133 if (acc != 0)
134 {
135 (*global_uo) << "measurement acc.=" << 2./acc << "%\n";
136 }
137 (*global_uo) << "TimeStep should be " << (TMGetTimer(3)+TMGetTimer(4)+TMGetTimer(5)+
138 TMGetTimer(6)+TMGetTimer(7)+TMGetTimer(8)+TMGetTimer(9)+TMGetTimer(11)+TMGetTimer(12)+TMGetTimer(20))/div*100 << "%.\n";
139
140 double tM=TMGetTimer(6);
141 double tK=TMGetTimer(5)-TMGetTimer(16);
142 double tC=TMGetTimer(8)+TMGetTimer(16);
143 double tS=TMGetTimer(14)+TMGetTimer(15);
144 double sum = tM+tK+tC+tS;
145 double sum2 = TMGetTimer(1);
146 double tR = sum2-sum;
147
148 if(sum2 != 0)
149 {
150 (*global_uo) << "CPU_time=" << TMGetTimer(1) << "\n";
151 (*global_uo) << "CPU%MKCSR=" << " " << tM/sum2*100 << "% " << tK/sum2*100 << "% " << tC/sum2*100
152 << "% " << tS/sum2*100 << "% " << tR/sum2*100 << "%\n";
153 }
154 /*
155 timings << TMGetTimer(1) << " " << tM/sum*100 << "% " << tK/sum*100 << "% " << tC/sum*100 << "% " << tS/sum*100 << "%\n" << flush;
156 timings2 << TMGetTimer(1) << " " << tM/sum2*100 << "% " << tK/sum2*100 << "% " << tC/sum2*100
157 << "% " << tS/sum2*100 << "% " << tR/sum2*100 << "%\n" << flush;
158 */
159 global_uo->ResetLocalMessageLevel();
160 #endif
161 }
162 //$ YV 2013-01-02: this is the implementation of timer functions for the kernel, the arrays can be directly accessed
163 #ifdef timeron
TMStartTimer(const int & i)164 inline void TMStartTimer(const int& i)
165 {
166 TMtstart(i) = GetClockTime();
167 }
168
TMStopTimer(const int & i)169 inline void TMStopTimer(const int& i)
170 {
171 TMtspent(i) += GetClockTime() - TMtstart(i);
172 }
173 #endif
TMGetTimer(int i)174 double TMGetTimer(int i)
175 {
176 if (TMtspent(i) < 1e-6) return 0;
177 else return TMtspent(i);
178 }
179
180 /*double Maximum(double a, double b, double c)
181 {
182 if (a>=b && a>=c) return a;
183 else if (b>=c) return b;
184 return c;
185 }
186 double Minimum(double a, double b, double c)
187 {
188 if (a<=b && a<=c) return a;
189 else if (b<=c) return b;
190 return c;
191 }*/
192
GetDiag(int n,double val)193 Matrix GetDiag(int n, double val)
194 {
195 int i;
196 Matrix m(n,n);
197 for (i=1; i<=n; i++)
198 m(i,i)=val;
199
200 return m;
201 }
202
IsEqual(const IVector & v1,const IVector & v2)203 int IsEqual(const IVector& v1, const IVector& v2)
204 {
205 if (v1.Length() != v2.Length()) return 0;
206
207 for (int i=1; i <= v1.Length(); i++)
208 if (v1(i) != v2(i)) return 0;
209
210 return 1;
211 }
212
Find(int i,const IVector & v)213 int Find(int i, const IVector& v)
214 {
215 for (int j=1; j <= v.Length(); j++)
216 {
217 if (v(j) == i) return j;
218 }
219 return 0;
220 }
221
222 //interpolate several values by means of cosine functions
InterpolateSmooth(double t,double t0,double t1,double A)223 double InterpolateSmooth(double t, double t0, double t1, double A)
224 {
225 double period = (t1-t0);
226 double tt = (t-t0)/period;
227 //double omega = MY_PI/period; //only half period of cosine fcn
228 //return A*(1.0-cos(omega*(t-t0)))/2.0;
229 return A*(10.*Cub(tt)-15.*To4(tt)+6*To5(tt));
230 }
231
232 //smooth interpolation over time t of function values f0,f1,f2,f3 at times t0,t1,t2,t3
InterpolateFnt2(double t,double t0,double t1,double t2,double t3,double f0,double f1,double f2,double f3)233 double InterpolateFnt2(double t, double t0, double t1, double t2, double t3, double f0, double f1, double f2, double f3)
234 {
235 if(t < t0) { return f0; }
236 else if(t < t1)
237 {
238 return InterpolateSmooth(t, t0, t1, f1-f0) + f0;
239 }
240 else if(t < t2)
241 {
242 return InterpolateSmooth(t, t1, t2, f2-f1) + f1;
243 }
244 else if(t < t3)
245 {
246 return InterpolateSmooth(t, t2, t3, f3-f2) + f2;
247 }
248 else return f3;
249 }
250
251 //smooth interpolation over time t of function values f0,f1,f2,f3 at times t0,t1,t2,t3
InterpolateFntArray(double t,TArray<double> t_vec,TArray<double> f_vec)252 double InterpolateFntArray(double t, TArray<double> t_vec, TArray<double> f_vec)
253 {
254 if (t_vec.Length() == 0 || t_vec.Length() != f_vec.Length()) return 0;
255
256 if(t < t_vec(1)) {return f_vec(1);}
257 for (int i=2; i <= t_vec.Length(); i++)
258 {
259 if (t < t_vec(i))
260 {
261 return InterpolateSmooth(t, t_vec(i-1), t_vec(i), f_vec(i)-f_vec(i-1)) + f_vec(i-1);
262 }
263 }
264 return f_vec.Last();
265 }
266
Mises(const Matrix & m)267 double Mises(const Matrix& m)
268 {
269 if (m.Getcols()!=3) return 0;
270
271 Matrix s = m-1./3.*m.Trace()*GetDiag(3);
272 return sqrt(3./2.*(s*s).Trace());
273
274 }
275
276 #ifdef gencnt
GenVecCnt()277 int GenVecCnt()
278 {
279 return *genvec;
280 }
GenMatCnt()281 int GenMatCnt()
282 {
283 return *genmat;
284 }
285 #endif
286
287
ApplyTransform(Vector & v,int mode,NumNLSys * nls)288 void SparseJacMat::ApplyTransform(Vector& v, int mode, NumNLSys* nls)
289 {
290 if (nls->TransformJacApply())
291 {
292 vtrans = v;
293 nls->ApplyTransform(vtrans,v,mode);
294 }
295 }
296
297 //********************** SPARSEJACMAT **********************************
Apply(const Vector & R,Vector & d,NumNLSys * nls)298 void SparseJacMat::Apply(const Vector& R, Vector& d, NumNLSys* nls)
299 {
300
301 TMStartTimer(15);
302
303 if (nls->UseSuperLU())
304 {
305
306 d = R;
307 //if (solvertype == PARDISOINVERSE)
308 //{
309 // pardisoinv->Apply(d);
310 //}
311 //else
312 //{
313 // superluinv->SuperLUApply(d);
314 //}
315 sparseinv->Apply(d);
316 }
317 else
318 {
319 //JG2013-10: not available open source
320 assert(0 && "Error: Band LU not available!");
321
322 //d.SetLen(J_vv.Getrows()+J_zz.Getrows());
323 ////mystatic Vector temp;
324 ////mystatic Vector temp2;
325 ////mystatic Vector Rsort;
326
327 //Rsort.SetLen(R.Length());
328 //if (resortsize)
329 //{
330 // //(*global_uo) << "resort R\n";
331 // for (int i = 1; i <= R.Length(); i++)
332 // Rsort(i) = R(resortvector(i));
333 //}
334 //else
335 //{
336 // Rsort = R;
337 //}
338
339
340 //d.SetAll(0);
341
342 //Vector R_v, R_z;
343 //R_v.LinkWith(Rsort,1,J_vv.Getcols());
344 //R_z.LinkWith(Rsort,1+J_vv.Getcols(),J_zz.Getcols());
345
346 //Vector q_v; q_v.LinkWith(d,1,J_vv.Getcols());
347 //Vector q_z; q_z.LinkWith(d,1+J_vv.Getcols(),J_zz.Getcols());
348
349 ////q_z = J_sch^(-1) * (R_z - (J_zv * (J_vv^(-1) * R_v)));
350
351 //temp = R_v;
352 //ApplyTransform(temp,0,nls);
353
354 //if (useband)
355 //{
356 // Band_LU_Bks0(J_vv_band,J_vv_band.Getrows(),lu,lu,J_vv_band_aux, J_vv_band_pivot, temp);
357 //}
358 //else
359 //{
360 // J_vv_band.LUBCKSUB(LUindx, temp);
361 //}
362
363 //ApplyTransform(temp,1,nls);
364
365
366 ////temp2 = R_v; //compare old and new Bks
367 ////Band_LU_Bks(J_vv_band,J_vv_band.Getrows(),lu,lu,J_vv_band_aux, J_vv_band_pivot, temp2);
368 ////(*global_uo) << "diff=" << (temp-temp2).MaxNorm() << "\n";
369
370 //Mult(J_zv,temp,temp2);
371 //temp = R_z;
372 //temp -= temp2;
373 //Mult(J_sch,temp,q_z);
374
375 ////q_v = J_vv^(-1) * (R_v - J_vz * q_z);
376
377 //Mult(J_vz,q_z,temp);
378 //temp2 = R_v;
379 //temp2 -= temp;
380
381 //q_v = temp2;
382 //ApplyTransform(q_v,0,nls);
383
384 //if (useband)
385 //{
386 // Band_LU_Bks0(J_vv_band,J_vv_band.Getrows(),lu,lu,J_vv_band_aux, J_vv_band_pivot, q_v);
387 //}
388 //else
389 //{
390 // J_vv_band.LUBCKSUB(LUindx, q_v);
391 //}
392
393 //ApplyTransform(q_v,1,nls);
394
395 //if (resortsize)
396 //{
397 // temp.SetLen(d.Length());
398 // temp = d;
399 // for (int i = 1; i <= d.Length(); i++)
400 // d(resortvector(i)) = temp(i);
401 //}
402 }
403 TMStopTimer(15);
404 }
405
406
407
Factorize(NumNLSys * nls)408 int SparseJacMat::Factorize(NumNLSys* nls)
409 {
410 //Daniel: comment out this function!!!
411 assert(0 && "Error: SparseJacMat::Factorize(NumNLSys* nls) not available any more!");
412
413 //double epszero = 1e-16; //optimize LUBCKSUB
414
415 ////(*global_uo) << "minv=" << minv_test << "\n";
416 ////mystatic Matrix temp1, temp2;
417 ////mystatic Vector tempcol;
418
419 //int n = J_vv.Getcols();
420 //int bw;
421 //int n_zz = J_zz.Getcols();
422
423 //if (LUcomputed) bw = LUcomputed_bw;
424 //else bw = J_vv.GetBandWidth();
425
426 ////(*global_uo) << "Jvv=" << J_vv << "\n";
427
428 //if (bw < (8*n)/10)
429 //{
430 // //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
431 // useband = 1;
432
433 // int col_band = bw*2-1;
434
435 // lu = bw-1;
436
437 // if (!(LUcomputed && nls->TransformJacApply()))
438 // {
439 // LUcomputed = 1;
440 // LUcomputed_bw = bw;
441 // TMStartTimer(27);
442
443 // J_vv_band_aux.SetSize(n, lu);
444
445 // J_vv.GetBandMatrix(J_vv_band, lu);
446 // //(*global_uo) << "Jzzband=" << J_vv_band << "\n";
447
448 // if (!n_vv_written) (*global_uo) << "LU-lalloc=" << (J_vv_band_aux.Getrows()*J_vv_band_aux.Getcols()+J_vv_band.Getrows()*J_vv_band.Getcols())*8./(1024.*1024.) << "MB\n";
449
450 // J_vv_band_pivot.SetLen(n);
451 // double signd;
452 // int rvb;
453 // rvb = Band_LU_Dec(J_vv_band,n,lu,lu,J_vv_band_aux, J_vv_band_pivot, signd);
454 // if (rvb == 0) {(*global_uo) << "banddec returned " << rvb << "\n"; return 0;}
455 // TMStopTimer(27);
456 // }
457
458 // //TMStartTimer(28);
459 // //J_sch = J_zz - J_zv * (J_vv^-1 * J_vz):
460 // mtemp1.SetSize(J_vv.Getrows(),J_vz.Getcols());
461 // tempcol.SetLen(J_vz.Getrows());
462 // for (int i = 1; i <= J_vz.Getcols(); i++)
463 // {
464 // J_vz.GetColVec(i,tempcol);
465 // if (!tempcol.IsZero(epszero))
466 // {
467 // ApplyTransform(tempcol,0,nls);
468 // Band_LU_Bks(J_vv_band,n,lu,lu,J_vv_band_aux, J_vv_band_pivot, tempcol);
469 // ApplyTransform(tempcol,1,nls);
470 // }
471 // mtemp1.SetColVec(tempcol,i);
472 // }
473 // //TMStopTimer(28);
474 //}
475 //else
476 //{
477 // //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
478 // //full version, J_vv_band is a full matrix and no banded matrix!!!
479 // useband = 0;
480
481 // TMStartTimer(27); //11.3% --> 14% with 500x500
482
483 // if (!(LUcomputed && nls->TransformJacApply()))
484 // {
485 // //(*global_uo) << "Jvv=" << J_vv << "\n";
486 // J_vv_band.CopyFrom(J_vv,1,1,n,n);
487
488 // //if (!n_vv_written) (*global_uo) << "J_vv=" << J_vv << "\n";
489 // //if (!n_vv_written) (*global_uo) << "J_vv_band=" << J_vv_band << "\n";
490
491 // LUcomputed = 1;
492 // LUcomputed_bw = bw;
493
494 // //double signd;
495 // int rvb = 0;
496 // rvb = J_vv_band.LU(LUindx, LUvv);
497 // if (rvb == 0) {(*global_uo) << "LU-decomposition returned " << rvb << "\n"; return 0;}
498 // }
499 // TMStopTimer(27);
500
501 // //J_sch = J_zz - J_zv * (J_vv^-1 * J_vz)
502 // //part: mtemp = (J_vv^-1 * J_vz):
503
504
505 // //TMStartTimer(28); //6% -> 3% with optimized LUBCKSUB
506 // mtemp1.SetSize(J_vv.Getrows(),J_vz.Getcols());
507 // tempcol.SetLen(J_vz.Getrows());
508 // //int cnt = 0;
509 // for (int i = 1; i <= J_vz.Getcols(); i++)
510 // {
511 // J_vz.GetColVec(i,tempcol);
512
513 // if (!tempcol.IsZero(epszero))
514 // {
515 // //cnt++;
516 // ApplyTransform(tempcol,0,nls);
517 // J_vv_band.LUBCKSUB(LUindx, tempcol);
518 // ApplyTransform(tempcol,1,nls);
519 // }
520 // mtemp1.SetColVec(tempcol,i);
521 // }
522 // //TMStopTimer(28);
523
524 // //(*global_uo) << "nonzero-cols=" << (double)cnt/((double)J_vz.Getcols()+1e-9)*100. << "% \n";
525
526 //}
527
528 //if (n_vv_written < 1) (*global_uo) << "n_vv=" << n << ", n_zz=" << J_zz.Getcols() << ", bw=" << bw << ", useband=" << useband << "\n";
529
530 ////J_sch = J_zz - J_zv * (J_vv^-1 * J_vz);
531 ////TMStartTimer(25); //dense: 12%, sparse version needs 0.5%
532 //mtemp2.SetSize(J_zz.Getrows(),J_zz.Getcols()); //needed???
533
534 //J_zvS.CopyFrom(J_zv); //sparse version
535 //Mult(J_zvS,mtemp1,mtemp2);
536
537 ////Mult(J_zv,mtemp1,mtemp2); //dense version
538
539 //J_sch = J_zz;
540 //J_sch -= mtemp2;
541 ////TMStopTimer(25);
542
543
544 //TMStartTimer(26); //<1%
545 //if (n_zz != 0)
546 //{
547 // if (!J_sch.Invert2()) {(*global_uo) << "Schur-complement could not be inverted!\n"; return 0;}
548 //}
549 //TMStopTimer(26);
550 //n_vv_written ++;
551
552 return 1;
553 }
554
555
556
557
558 //********************** NUMNLSYS**********************************
UOfull(int message_level,int output_prec) const559 UserOutput& NumNLSys::UOfull(int message_level, int output_prec) const
560 {
561 uo.SetLocalMessageLevel(message_level);
562 uo.SetOutputPrec(output_prec);
563 return uo;
564 };
565
SetSolver(NumNLSolver * s)566 void NumNLSys::SetSolver(NumNLSolver* s)
567 {
568 solver = s;
569 }
570
NLS_GetJacCol() const571 int NumNLSys::NLS_GetJacCol() const
572 {
573 if (solver->SymmetricJacobian())
574 {
575 return jaccol;
576 }
577 else {return 0;}
578 }
579
UseSparseSolver() const580 int NumNLSys::UseSparseSolver() const
581 {
582 //if (TransformJacApply()) return 1;
583 //return 0;
584 return solver->UseSparseSolver();
585 }
586
SolveUndeterminedSystem() const587 int NumNLSys::SolveUndeterminedSystem() const
588 {
589 return solver->SolveUndeterminedSystem();
590 }
591
SolveUndeterminedSystem()592 int& NumNLSys::SolveUndeterminedSystem()
593 {
594 return solver->SolveUndeterminedSystem();
595 }
596
EstimatedConditionNumber() const597 double NumNLSys::EstimatedConditionNumber() const
598 {
599 return solver->EstimatedConditionNumber();
600 }
601
EstimatedConditionNumber()602 double& NumNLSys::EstimatedConditionNumber()
603 {
604 return solver->EstimatedConditionNumber();
605 }
606
UseSparseJac() const607 int NumNLSys::UseSparseJac() const
608 {
609 //return 0;
610 return solver->UseSparseSolver();
611 }
612
UseSuperLU() const613 int NumNLSys::UseSuperLU() const
614 {
615 return 1;
616 }
617
Jacobian(Matrix & m,Vector & x)618 void NumNLSys::Jacobian(Matrix& m, Vector& x)
619 {
620 TMStartTimer(3);
621 int i,j;
622 //mystatic Vector f1;
623 //mystatic Vector f2;
624 f1.SetLen(x.GetLen());
625 f2.SetLen(x.GetLen());
626
627 m.SetSize(x.GetLen(),x.GetLen());
628
629 double numdiffepsi = solver->NumDiffepsi();
630 NLS_SetJacCol(0);
631 NLS_SetJacFullNewton(1);
632 fullnewtoncnt++;
633
634 if (!solver->SymmetricJacobian())
635 {
636 NLF(x,f2);
637 }
638 double eps;
639 double xstore;
640 (*global_uo) << "len=" << m.Getcols() << "\n";
641 for (i = 1; i <= m.Getcols(); i++)
642 {
643 eps = numdiffepsi*Maximum(1e-2,fabs(x(i)));
644 if (solver->SymmetricJacobian())
645 {
646 //if (TIstages==1)
647 //{
648 NLS_SetJacCol(i);
649 //}
650
651 xstore = x(i);
652 x(i) += eps;
653 NLF(x,f1);
654 //(*global_uo) << "f1=" << f1 << "\n";
655
656 x(i) -= 2*eps;
657 NLF(x,f2);
658
659 //(*global_uo) << "f2=" << f2 << "\n";
660 x(i) = xstore;
661 }
662 else
663 {
664 xstore = x(i);
665 x(i) += 2*eps;
666 NLF(x,f1);
667 x(i) = xstore;
668 }
669
670 for (j=1; j<=f1.GetLen(); j++)
671 {
672 m(j,i)=0.5/eps*(f1(j)-f2(j));
673 }
674 }
675
676 NLS_SetJacFullNewton(0);
677 NLS_SetJacCol(0);
678 TMStopTimer(3);
679 }
680
681 int lastmatrix = 0;
682 Matrix lastm;
683
Factorize(Matrix & minv,SparseMatrix & sminv,SparseJacMat & sparseminv)684 int NumNLSolver::Factorize(Matrix& minv, SparseMatrix& sminv, SparseJacMat& sparseminv)
685 {
686 //Matrix mm = sminv.GetMatrix();
687 //PrintMatrix01(mm);
688 //global_uo.InstantMessageText("NumSolver::Factorize begin");
689
690 global_uo->SaveLocalMessageLevel();
691 global_uo->SetLocalMessageLevel(UO_LVL_dbg1);
692
693 TMStartTimer(14);
694 int rv = 0;
695 if (nls->UseSparseSolver())
696 {
697 if (nls->UseSparseJac()) //use sminv
698 {
699 if (nls->UseSuperLU())
700 {
701 double facttime;
702 delete sparseminv.sparseinv;
703 sparseminv.sparseinv = new SparseInverse(sminv);
704 sparseminv.sparseinv->SetSparseMatrix(sminv);
705 rv = 1;
706 if (sminv.Getcols()>10000 && global_uo->PrintMsg()) (*global_uo) << "Factorize Matrix with sparse solver .. ";
707 facttime = -GetClockTime();
708 rv = sparseminv.sparseinv->Factorize();
709 facttime += GetClockTime();
710 if (sminv.Getcols()>10000 && global_uo->PrintMsg()) (*global_uo) << "done, factorization time " << facttime << "s\n";
711 sparseminv.sparseinv->IsFirstStep()=0;
712
713 }
714 else
715 {
716 //JG2013-10: functions not available open source
717 assert(0 && "Error: Band Sparse Factorize not available any more!");
718 // //(*global_uo) << "factorize changed\n";
719 // sparseminv.resortsize = nls->GetResortSize();
720 // int rs = nls->GetResortVector().Length();
721 // sparseminv.resortvector.SetLen(rs);
722
723 // for (int i=1; i <= rs; i++)
724 // {
725 // sparseminv.resortvector(i) = nls->GetResortVector()(i);
726 // }
727 // /*
728 // (*global_uo) << "J=" << sminv << "\n";
729 // (*global_uo) << "resortsize=" << sparseminv.resortsize << "\n";
730 // (*global_uo) << "rs=" << rs << "\n";
731 // (*global_uo) << "rsvector=" << sparseminv.resortvector << "\n";
732 // */
733 // if (nls->GetResortSize())
734 // {
735 // sparseminv.J_vv.CopyFrom(sminv,1,1,bandsize,bandsize,nls->GetResortVector());
736 // sparseminv.J_vz.CopyFrom(sminv,1,1+bandsize,bandsize,sminv.Getcols(),nls->GetResortVector());
737 // sparseminv.J_zv.CopyFrom(sminv,1+bandsize,1,sminv.Getrows(),bandsize,nls->GetResortVector());
738 // sparseminv.J_zz.CopyFrom(sminv,1+bandsize,1+bandsize,
739 // sminv.Getrows(),sminv.Getcols(),nls->GetResortVector());
740 // }
741 // else
742 // {
743 // sparseminv.J_vv.CopyFrom(sminv,1,1,bandsize,bandsize);
744 // //global_uo.InstantMessageText("NumSolver::Factorize begin2");
745 // sparseminv.J_vz.CopyFrom(sminv,1,1+bandsize,bandsize,sminv.Getcols());
746 // sparseminv.J_zv.CopyFrom(sminv,1+bandsize,1,sminv.Getrows(),bandsize);
747 // sparseminv.J_zz.CopyFrom(sminv,1+bandsize,1+bandsize,
748 // sminv.Getrows(),sminv.Getcols());
749
750 // //(*global_uo) << "sminv=" << sminv << "\n";
751 // //(*global_uo) << "J_vv=" << sparseminv.J_vv << "\n";
752 // //(*global_uo) << "J_vz=" << sparseminv.J_vz << "\n";
753 // //(*global_uo) << "J_zv=" << sparseminv.J_zv << "\n";
754 // //(*global_uo) << "J_zz=" << sparseminv.J_zz << "\n";
755
756
757 // }
758
759 // rv = sparseminv.Factorize(nls);
760 }
761 }
762 else
763 {
764 //JG2013-10: case cannot happen
765 assert(0 && "Error: Band Sparse Factorize not available any more!");
766 //// bandsize = minv.Getcols(); //for testing
767 //sparseminv.resortsize = nls->GetResortSize();
768 //if (nls->GetResortSize())
769 //{
770 // //TMStartTimer(24); //1% time lost
771
772 // int rs = nls->GetResortVector().Length();
773 // sparseminv.resortvector.SetLen(rs);
774 // for (int i=1; i <= rs; i++)
775 // {
776 // sparseminv.resortvector(i) = nls->GetResortVector()(i);
777 // }
778
779
780 // //(*global_uo) << "sm-diff=" << (sm.GetMatrix()-minv).MaxNorm() << "\n";
781
782 // //original:
783 // //SparseMatrix sm(minv);
784 // //sparseminv.J_vv.CopyFrom(minv,1,1,bandsize,bandsize,nls->GetResortVector());
785 //
786 // sparseminv.mtemp1.CopyFrom(minv,1,1,bandsize,bandsize,nls->GetResortVector()); //1% time lost
787 // sparseminv.J_vv.CopyFrom(sparseminv.mtemp1);
788
789 // //sparseminv.J_vv.CopyFrom(minv,1,1,bandsize,bandsize,nls->GetResortVector());
790 // sparseminv.J_vz.CopyFrom(minv,1,1+bandsize,bandsize,minv.Getcols(),nls->GetResortVector());
791 // sparseminv.J_zv.CopyFrom(minv,1+bandsize,1,minv.Getrows(),bandsize,nls->GetResortVector());
792 // sparseminv.J_zz.CopyFrom(minv,1+bandsize,1+bandsize,minv.Getrows(),minv.Getcols(),nls->GetResortVector());
793
794 // //TMStopTimer(24);
795
796 //}
797 //else
798 //{
799 // TMStartTimer(24);
800
801 // sparseminv.J_vv.CopyFrom(minv,1,1,bandsize,bandsize);
802 // sparseminv.J_vz.CopyFrom(minv,1,1+bandsize,bandsize,minv.Getcols());
803 // sparseminv.J_zv.CopyFrom(minv,1+bandsize,1,minv.Getrows(),bandsize);
804 // sparseminv.J_zz.CopyFrom(minv,1+bandsize,1+bandsize,
805 // minv.Getrows(),minv.Getcols());
806
807 // TMStopTimer(24);
808
809 //}
810 ////sparseminv.minv_test = minv;
811 ////sparseminv.minv_test.Invert2();
812
813 //rv = sparseminv.Factorize(nls);
814 }
815 }
816 else
817 {
818 if (nls->UseSparseJac()) //use sminv
819 {
820 minv.CopyFrom(sminv,1,1,sminv.Getrows(),sminv.Getcols());
821 }
822
823 rv = minv.Invert2();
824
825 }
826 TMStopTimer(14);
827
828 global_uo->ResetLocalMessageLevel();
829
830 return rv;
831 }
832
833
834 //********************** APPLYJAC **********************************
ApplyJac(const Vector & f,Vector & res)835 int NumNLSolver::ApplyJac(const Vector& f, Vector& res)
836 {
837 SaveJac* jac;
838 int chj = ChooseJac();
839 //(*global_uo) << "nlsinfo=" << nonlinsolveinfo
840 // << ", nls1=" << sjac[0].nlsinfo << ", nls2=" << sjac[1].nlsinfo << "\n";
841 if (chj != 0)
842 {
843 jac = &(sjac[chj-1]);
844 Mult(jac->oldjacmat,f,res);
845 return 1;
846 }
847 else return 0;
848 }
849
850
851 //********************** NUMNLSOLVER**********************************
NLSolve(Vector & x0)852 int NumNLSolver :: NLSolve(Vector& x0)
853 {
854 if (SolverPrintsDetailedOutput())
855 {
856 SolverUO() << "-->";
857 if (ModifiedNewton())
858 {
859 SolverUO() << "Modified Newton:";
860 }
861 else
862 {
863 SolverUO() << "Full Newton:";
864 }
865 }
866
867 TMStartTimer(10);
868 //(*global_uo) << "NLSolve\n";
869 //(*global_uo) << "genvec=" << genvec << ", ";
870 //(*global_uo) << "genmat=" << genmat << "\n";
871 //(*global_uo) << "nls-inf=" << this->nonlinsolveinfo << "\n";
872 //(*global_uo) << " xstart0=" << x0 << "\n";
873
874 jaccount = 0;
875 nls->SetSolver(this);
876 //mystatic Vector x0start;
877 x0start = x0;
878 nls->NLS_SetJacCol(0);
879
880 //mystatic Vector f;
881 f.SetLen(x0.GetLen());
882 //mystatic Vector xd;
883 xd.SetLen(x0.GetLen());
884
885 const double MAXERR = 1e100;
886 const double MAXERRINC = 1e10; //$JG 2012-11-07 changed from 1e10 to 1e25 because of problems with Cable2D element in first step
887 int it = 0;
888
889 nls->NLF(x0,f);
890
891
892 double errinit = f.MaxNorm();
893 double err = errinit;
894 if (errinit < absoluteaccuracy)
895 {
896 errinit = absoluteaccuracy;
897 }
898 double error_goal = relativeaccuracy*errinit; // error_goal = relativeaccuracy * min(absoluteaccuracy, errinit)
899
900 if (SolverPrintsDetailedOutput())
901 {
902 mystr str(" goal=");
903 str += mystr(error_goal);
904 str += mystr(" (err=");
905 str += mystr(err);
906 str += mystr(", rel_acc=");
907 str += mystr(relativeaccuracy);
908 str += mystr(", abs_acc=");
909 str += mystr(absoluteaccuracy);
910 str += mystr(", [goal=max(err,abs_acc)rel_acc, err=MaxNorm(residual)])\n");
911
912 SolverUO() << str;
913
914 if (nls->GetOptions()->LoggingOptions()->SolverNewtonIterationResidualVector())
915 {
916 SolverUO() << "residual vector = " << f << "\n";
917 }
918 }
919
920 double lasterr = err;
921 int jacsingular = 0;
922 int stopnewton = 0;
923
924 double maxcontractivity = 0;
925 contractivity = 0.5;
926 Matrix* minv; //full matrix
927 SparseMatrix* sminv; //sparse matrix build
928 SparseJacMat* sparseminv; //sparse matrix factorize
929
930 //initialize for full newton:
931 SaveJac* jac = &sjac[0];
932 minv = &(jac->oldjacmat);
933 sminv = &(jac->oldsjacmat);
934 sparseminv = &(jac->sparsejacmat);
935
936 if (ModifiedNewton())
937 {
938 int chj = ChooseJac();
939 //chj=0; //every timestep a Jacobian
940 if (chj != 0)
941 {
942 jac = &(sjac[chj-1]);
943 }
944 else
945 {
946 if (sjac[0].oldjac == 0) {jac = &sjac[0]; }
947 else if (sjac[1].oldjac == 0) {jac = &sjac[1];}
948 else if ((fabs(sjac[0].nlsinfo) >= fabs(nonlinsolveinfo) || fabs(sjac[1].nlsinfo) >= fabs(nonlinsolveinfo))
949 && sjac[0].lastbuilt) {jac = &sjac[1];}
950 else if ((fabs(sjac[0].nlsinfo) >= fabs(nonlinsolveinfo) || fabs(sjac[1].nlsinfo) >= fabs(nonlinsolveinfo))
951 && sjac[1].lastbuilt) {jac = &sjac[0];}
952 //else if ((sjac[1].oldjacage < sjac[0].oldjacage && sjac[1].lastbuilt)) {jac = &sjac[0];}
953 //else if ((sjac[0].oldjacage < sjac[1].oldjacage && sjac[0].lastbuilt)) {jac = &sjac[1];}
954 else if ((fabs(sjac[1].nlsinfo - nonlinsolveinfo) > fabs(sjac[0].nlsinfo - nonlinsolveinfo))) {jac = &sjac[1]; }
955 else {jac = &sjac[0];}
956
957 jac->oldjac = 0;
958 }
959
960 minv = &(jac->oldjacmat);
961 sminv = &(jac->oldsjacmat);
962 sparseminv = &(jac->sparsejacmat);
963
964 if (jac->oldjac == 0 || jac->oldjacage > jac->maxjacage || jac->oldjacsize != x0.GetLen())
965 {
966
967 // report on the reason to update jacobian
968 if (SolverPrintsDetailedOutput())
969 {
970 mystr jacobian_update_info_sparsity = (nls->UseSparseJac() ? mystr("Sparse ") : mystr("Full "));
971 mystr jacobian_update_info = jacobian_update_info_sparsity + mystr("Jacobian (#") + mystr(jaccount) + mystr(") is computed");
972 if (jac->oldjac == 0)
973 {
974 jacobian_update_info += mystr(".\n");
975 }
976 else if (jac->oldjacage > jac->maxjacage)
977 {
978 jacobian_update_info += mystr(", since the age of the Jacobian was greater than ") + mystr(jac->maxjacage) + mystr(".\n");
979 }
980 else //jac->oldjacsize != x0.GetLen()
981 {
982 jacobian_update_info += mystr(", since the length of the solution vector did not fit the size of the Jacobian.\n");
983 }
984 SolverUO() << jacobian_update_info;
985 }
986
987 AssembleJacobian(minv, sminv, x0);
988
989 /*
990 nls->Jacobian(*sminv,x0);
991 nls->Jacobian(*minv,x0);
992 //(*global_uo) << "Jac=" << *minv;
993 //(*global_uo) << "Jac_s=" << sminv->GetMatrix();
994 //(*global_uo) << "Jac-diff=" << *minv-sminv->GetMatrix();
995 (*global_uo) << "N-Jac-diff=" << (*minv-sminv->GetMatrix()).MaxNorm() << "\n";
996 */
997
998 if (!Factorize(*minv,*sminv,*sparseminv))
999 {
1000 jacsingular = 1000;
1001 };
1002
1003 jac->oldjac = 1;
1004 jac->nlsinfo = nonlinsolveinfo;
1005 jac->oldjacage = 0;
1006 jac->oldjacsize = x0.GetLen();
1007 sjac[0].lastbuilt = 0;
1008 sjac[1].lastbuilt = 0;
1009 jac->lastbuilt = 1;
1010 }
1011 else
1012 {
1013 jac->oldjacage++;
1014 }
1015
1016 }
1017
1018 int n_modstopped = 0;
1019 stopmnr = 0;
1020 int contit = 0;
1021
1022 while ((it < MaxNewtonSteps() && err >= error_goal && jacsingular < 2 && !stopnewton) || it==0)
1023 {
1024 it++;
1025 contit++;
1026 int maxmodnewtonsteps1 = MaxModNewtonSteps();
1027 int maxmodnewtonsteps2 = MaxModNewtonSteps()+MaxRestartNewtonSteps();
1028
1029 double low_contractivity_tolerance = 0.7;
1030 double high_contractivity_tolerance = 2.;
1031 if (((it == maxmodnewtonsteps1) || (contit>2 && contractivity > low_contractivity_tolerance)) && !stopmnr && ModifiedNewton()) //contractivity > 0.7 is optimum for many cases
1032 {
1033 // (*global_uo) << "maxmodnewtonsteps1=" << maxmodnewtonsteps1 << ",contr=" << contractivity << "\n";
1034 double memo_contractivity = contractivity;
1035 if (n_modstopped==0)
1036 {
1037 x0=x0start;
1038 contractivity = 0.5;
1039 }
1040 if ((contractivity >= high_contractivity_tolerance && n_modstopped>0) || n_modstopped>2)
1041 {
1042 if (SolverPrintsDetailedOutput())
1043 {
1044 SolverUO() << mystr("Newton error was set to MAXERR=")+mystr(MAXERR)+mystr(", since ");
1045 if (contractivity >= high_contractivity_tolerance)
1046 {
1047 SolverUO() << mystr("contractivity=")+mystr(contractivity)+mystr(" exceeded highly critical value of ")+mystr(high_contractivity_tolerance)+mystr(".\n");
1048 }
1049 else
1050 {
1051 SolverUO() << mystr("Jacobian has been updated more than twice, but still ");
1052 if (contractivity > low_contractivity_tolerance)
1053 {
1054 SolverUO() << mystr("contractivity=")+mystr(memo_contractivity)+mystr(" exceeded critical value of ")+mystr(low_contractivity_tolerance)+mystr(".\n");
1055 }
1056 else
1057 {
1058 SolverUO() << mystr("iterations=")+mystr(it)+mystr(" exceeded critical value of ")+mystr(maxmodnewtonsteps1)+mystr(".\n");
1059 }
1060 }
1061 }
1062 err=MAXERR;
1063 }
1064 else
1065 {
1066 contit = 0;
1067 n_modstopped++;
1068
1069
1070 // (*global_uo) << "jac1\n";
1071
1072 if (SolverPrintsDetailedOutput())
1073 {
1074 mystr jacobian_update_info_sparsity = (nls->UseSparseJac() ? mystr("Sparse ") : mystr("Full "));
1075 mystr jacobian_update_info = jacobian_update_info_sparsity + mystr("Jacobian (# ") + mystr(jaccount) + mystr(") is computed");
1076 if (it == maxmodnewtonsteps1)
1077 {
1078 jacobian_update_info += mystr(", since iterations=") + mystr(it) + mystr(" exceeded critical value of ") + mystr(maxmodnewtonsteps1) + mystr(".\n");
1079 }
1080 else
1081 {
1082 jacobian_update_info += mystr(", since contractivity ") + mystr(memo_contractivity) + mystr(" exceeded critical value of ") + mystr(low_contractivity_tolerance) + mystr(", and contit is greater than 2.\n");
1083 }
1084 SolverUO() << jacobian_update_info;
1085 }
1086
1087 AssembleJacobian(minv, sminv, x0);
1088
1089 nls->NLF(x0,f);
1090 lasterr = f.MaxNorm();
1091
1092 if (!Factorize(*minv,*sminv,*sparseminv))
1093 {
1094 jacsingular = 1000;
1095 };
1096
1097
1098 jac->oldjac = 1;
1099 jac->nlsinfo = nonlinsolveinfo;
1100 jac->oldjacage = 0;
1101 jac->oldjacsize = x0.GetLen();
1102 sjac[0].lastbuilt = 0;
1103 sjac[1].lastbuilt = 0;
1104 jac->lastbuilt = 1;
1105 }
1106 }
1107 if (ModifiedNewton() && (it > maxmodnewtonsteps2 || err == MAXERR) && !stopmnr)
1108 {
1109 if (SolverPrintsDetailedOutput())
1110 {
1111 SolverUO() << mystr("Switch from Modified to Full Newton, since ");
1112 if (it > maxmodnewtonsteps2)
1113 {
1114 SolverUO() << mystr("iterations=")+mystr(it)+mystr(" exceed highly critical value of ")+mystr(maxmodnewtonsteps2)+mystr(".\n");
1115 }
1116 else
1117 {
1118 SolverUO() << mystr("Newton error=MAXERR=")+mystr(MAXERR)+mystr(".\n");
1119 }
1120 }
1121
1122 stopmnr = 1;
1123 x0 = x0start;
1124 nls->NLF(x0,f);
1125 jac->oldjac = 0;
1126 jac->jacfailedcnt++;
1127 }
1128 //if (jac->jacfailedcnt > 20) {jac->maxjacage = jac->maxjacage/2+1; jac->jacfailedcnt = 0;}
1129
1130 if (!ModifiedNewton() || stopmnr)
1131 {
1132 AssembleJacobian(minv, sminv, x0);
1133
1134 // factorize Jacobian & solve system
1135 contractivity = 0;
1136 fulljaccnt++;
1137 //mystatic IVector iv;
1138 int rv;
1139
1140 if (nls->UseSparseSolver())
1141 {
1142 rv = Factorize(*minv,*sminv,*sparseminv);
1143 if (rv) sparseminv->Apply(f,xd,nls);
1144 }
1145 else
1146 {
1147 if (nls->SolveUndeterminedSystem())
1148 {
1149 int minv_rank = 0;
1150 double rcond = 1./nls->EstimatedConditionNumber();
1151 rv = minv->UndeterminedSystemSolve(f, xd, rcond, minv_rank);
1152 int minv_dim = min(minv->Getrows(), minv->Getcols());
1153 if (SolverPrintsDetailedOutput() && minv_rank < minv_dim)
1154 {
1155 SolverUO() << "WARNING: Estimated rank of Jacobian (" << minv_rank << ") smaller than dimension (" << minv_dim << ")\n";
1156 }
1157 }
1158 else
1159 {
1160 rv = minv->Solve(f,xd);
1161 }
1162 }
1163
1164
1165 if (0) // old branch
1166 {
1167 rv = minv->LU(iv);
1168 if (rv < 0) rv=0;
1169 xd=f;
1170 if (rv) minv->LUBCKSUB(iv,xd);
1171 }
1172
1173 if (!rv)
1174 {
1175 x0 = x0start;
1176 xd.SetAll(0);
1177 nls->NLF(x0,f);
1178 jacsingular+=1;
1179 };
1180 }
1181 else
1182 {
1183 //TMStartTimer(24);
1184 if (nls->UseSparseSolver())
1185 {
1186 sparseminv->Apply(f,xd,nls);
1187 }
1188 else
1189 {
1190 TMStartTimer(15);
1191 Mult(*minv,f,xd);
1192 TMStopTimer(15);
1193 //(*global_uo) << "xdf=" << xd << "\n";
1194 }
1195
1196 //TMStopTimer(24);
1197 }
1198
1199 double trust = 1;
1200
1201 if (trustregion)
1202 {
1203 Vector x;
1204 double lasterr, acterr;
1205 trust = trustregiondiv;
1206 lasterr = 1e50;
1207
1208 double ftrust = 0;
1209
1210 while (trust <= 2)
1211 {
1212 x=x0-trust*xd; //Vector generation!!!!
1213 nls->NLF(x,f);
1214
1215 acterr = f.MaxNorm();
1216
1217 if (acterr < lasterr) {ftrust = trust; lasterr = acterr;}
1218
1219 trust += trustregiondiv;
1220 }
1221
1222 trust = ftrust;
1223 if (ftrust == 0) {trust = 1;}
1224
1225 }
1226
1227 //trust = 0.8;//0.8 work best, DASSL uses 0.75 - but this leads to many Jacobians sometimes!!!; 1-jac->oldjacage*0.025 ??;
1228 if (trust != 1) {xd.Mult(trust);}
1229
1230 x0 -= xd;
1231 nls->NLF(x0,f);
1232 err = f.MaxNorm();
1233
1234 if (SolverPrintsDetailedOutput())
1235 {
1236 mystr str("Nit=");
1237 str += mystr(it);
1238 str += mystr(": err=");
1239 str += mystr(err);
1240 if (ModifiedNewton() && !stopmnr)
1241 {
1242 str += mystr(", contractivity=");
1243 str += mystr(contractivity);
1244 }
1245 str += mystr("\n");
1246 SolverUO() << str;
1247
1248 if (nls->GetOptions()->LoggingOptions()->SolverNewtonIterationSolutionVector())
1249 {
1250 SolverUO() << "solution vector = " << xd << "\n";
1251 }
1252 if (nls->GetOptions()->LoggingOptions()->SolverNewtonIterationResidualVector())
1253 {
1254 SolverUO() << "residual vector = " << f << "\n";
1255 }
1256 //EK & PG 2012-05-18 at the end of the computation no recomputation of jacobian necessary - see comment below
1257 //if (nls->GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiMatrix())
1258 //{
1259 // if (nls->UseSparseJac())
1260 // {
1261 // //SparseMatrix m; nls->Jacobian(m,f); //PG 2012-05-18: recalculation of sparse jacobian lead to different results for (SolverNewtonIterationJacobiMatrix()==0/1)
1262 // Matrix m = sminv->GetMatrix();
1263 // SolverUO() << "sparse jacobi matrix = " << m << "\n";
1264 // }
1265 // else
1266 // {
1267 // //Matrix m; nls->Jacobian(m,f); //PG 2012-05-18: recalculation not necessary, however no difference behaviour (as in sparse case) was encountered
1268 // Matrix m(*minv);
1269 // SolverUO() << "jacobi matrix = " << m << "\n";
1270 // }
1271 //}
1272 }
1273
1274 if (lasterr != 0)
1275 {
1276 if (contractivity == 0)
1277 {
1278 contractivity = err/lasterr;
1279 }
1280 else
1281 {
1282 contractivity = sqrt(contractivity*err/lasterr);
1283 }
1284 }
1285 else
1286 {
1287 contractivity = 0.1;
1288 }
1289
1290 if (it>4) maxcontractivity=Maximum(maxcontractivity,contractivity);
1291 lasterr = err;
1292
1293 //MaSch 2013-02-04: temporarily changed the MAXERRINC-test to avoid comparison to absolute numbers in the first Newton step (where x0start = 0)
1294 //due to problems with very small and thin (length ~ 1e-6 m) ANCFCable2D elements
1295 //(test case: finite bending deformation by constant load in x-direction with clamped end, r' in x-dir gets very large in the first step (scaling linearly with force/density),
1296 //such that x0.MaxNorm() > MAXERRINC*(x0start.MaxNorm()+1.) = MAXERRINC in the first Newton step
1297 //static computation without problems;
1298
1299 //if (!x0.IsValid(MAXERR) || IsNaN(err) || (err >= MAXERR) || (err <= -MAXERR) || x0.MaxNorm() > MAXERRINC*(x0start.MaxNorm()+1.)) //err is Nan or out of range
1300 if (!x0.IsValid(MAXERR) || IsNaN(err) || (err >= MAXERR) || (err <= -MAXERR) || (x0start.MaxNorm()==0. ? false : (x0.MaxNorm() > MAXERRINC*(x0start.MaxNorm()))) ) //err is Nan or out of range
1301 {
1302 if(SolverPrintsDetailedOutput())
1303 {
1304 mystr str("err was set to MAXERR=");
1305 str+=mystr(MAXERR);
1306 str+=mystr(", since\n");
1307 if (!x0.IsValid(MAXERR)) str+=mystr(" * max-norm of x0 exceeds critical value of MAXERR=")+mystr(MAXERR)+mystr("\n");
1308 if (IsNaN(err)) str+=mystr(" * max-norm of residual is nan\n");
1309 if (err >= MAXERR || err <= -MAXERR) str+=mystr(" * max-norm of residual exceeds critical value of MAXERR=")+mystr(MAXERR)+mystr("\n");
1310 if (x0.MaxNorm() > MAXERRINC*(x0start.MaxNorm()+1.)) str+=mystr(" * max-norm of x0 exceeds critical value of MAXERRINC*(x0start.MaxNorm()+1.)=")+mystr(MAXERRINC*(x0start.MaxNorm()+1.))+mystr("\n");
1311 SolverUO() << str;
1312 }
1313 err = MAXERR;
1314
1315 if (!ModifiedNewton() || stopmnr)
1316 {
1317 stopnewton = 1;
1318
1319 if(SolverPrintsDetailedOutput())
1320 {
1321 mystr str("Full Newton method stopped!\n");
1322 SolverUO() << str;
1323 }
1324 }
1325 else
1326 {
1327 }
1328 }
1329 }
1330
1331 newtonits = it;
1332 contractivity = maxcontractivity;
1333 TMStopTimer(10);
1334
1335 if (!x0.IsValid(MAXERR)) //values of solution are out of range
1336 {
1337 error_msg += "Solution diverged! ";
1338 if (SolverPrintsDetailedOutput())
1339 {
1340 SolverUO() << " error message=Solution diverged! " << "\n";
1341 }
1342 if (jacsingular != 0)
1343 {
1344 error_msg += "Jacobian is singular! ";
1345 if (SolverPrintsDetailedOutput())
1346 {
1347 SolverUO() << " error message=Jacobian is singular! " << "\n";
1348 }
1349 }
1350 return 0;
1351 }
1352
1353 if (err == MAXERR) //err is NaN
1354 {
1355 error_msg += "Residual went to infinity! ";
1356 if (SolverPrintsDetailedOutput()) {SolverUO() << " error message=Residual went to infinity! " << "\n";}
1357 if (jacsingular != 0)
1358 {
1359 error_msg += "Jacobian is singular! ";
1360 if (SolverPrintsDetailedOutput()) {SolverUO() << " error message=Jacobian is singular! " << "\n";}
1361 }
1362 x0 = x0start;
1363 return 0;
1364 }
1365
1366 if ((err < error_goal) && !jacsingular && !stopnewton)
1367 {
1368 return 1;
1369 }
1370 else
1371 {
1372 if (jacsingular != 0)
1373 {
1374 error_msg += "Jacobian is singular! ";
1375 if (SolverPrintsDetailedOutput()) {SolverUO() << " error message=Jacobian is singular! " << "\n";}
1376 }
1377 else if (err >= error_goal)
1378 {
1379 error_msg += "Error tolerance not reached (with Full Newton)!\n";
1380 error_msg += mystr("Error = ")+mystr(err)+mystr(", err-goal=") + mystr(error_goal);
1381 if (SolverPrintsDetailedOutput()) {SolverUO() << " error message=Error tolerance not reached (with Full Newton)!" << "\n";}
1382 }
1383 else
1384 {
1385 error_msg += "Newton method not successful! ";
1386 if (SolverPrintsDetailedOutput()) {SolverUO() << " error message=Newton method not successful! " << "\n";}
1387 }
1388
1389 x0 = x0start;
1390 return 0;
1391 }
1392
1393 }
1394
1395 //$ PG 2013-11-6: assemble Jacobian (and print some information), set jaccount++
AssembleJacobian(Matrix * minv,SparseMatrix * sminv,Vector & x0)1396 void NumNLSolver::AssembleJacobian(Matrix* minv, SparseMatrix* sminv, Vector& x0)
1397 {
1398 if (nls->UseSparseJac())
1399 {
1400 nls->Jacobian(*sminv,x0);
1401 jaccount++;
1402 if (GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiMatrix() && SolverPrintsDetailedOutput())
1403 {
1404 OutputJacobian(*sminv);
1405 }
1406 if (GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiCondition() && SolverPrintsDetailedOutput())
1407 {
1408 OutputConditionNumber(*sminv);
1409 }
1410 }
1411 else
1412 {
1413 nls->Jacobian(*minv,x0);
1414 jaccount++;
1415 if (GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiMatrix() && SolverPrintsDetailedOutput())
1416 {
1417 OutputJacobian(*minv);
1418 }
1419 if (GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiCondition() && SolverPrintsDetailedOutput())
1420 {
1421 OutputConditionNumber(*minv);
1422 }
1423 }
1424 }
1425
1426 //$ PG 2013-10-17: calculate condition number of jacobi matrix, IS SLOW, should only be used for designing elements, constraints, for adjusting models, or for general debugging purpose
OutputConditionNumber(const Matrix & M) const1427 void NumNLSolver::OutputConditionNumber(const Matrix& M) const
1428 {
1429 double rcond;
1430 int info = M.EstimateReciprocalConditionNumber(rcond);
1431
1432 if (info < 0)
1433 {
1434 SolverUO() << "ERROR: Illegal value in matrix at position " << -info << " of the array!\n";
1435 }
1436 else if (info > 0)
1437 {
1438 SolverUO() << "ERROR: Zero division in LU factorization. M(" << info << "," << info << ") = 0, matrix is singular!\n";
1439 }
1440 else //if (info==0)
1441 {
1442 SolverUO() << "kappa = " << 1./rcond << "\n";
1443 }
1444 }
1445
1446 //$ PG 2013-10-17: calculate condition number of jacobi matrix, IS SLOW, should only be used for designing elements, constraints, for adjusting models, or for general debugging purpose
OutputConditionNumber(const SparseMatrix & M) const1447 void NumNLSolver::OutputConditionNumber(const SparseMatrix& M) const
1448 {
1449 SolverUO() << "WARNING: Calculation of condition number for sparse matrices uses routine for full matrices!\n";
1450 OutputConditionNumber(M.GetMatrix());
1451 }
1452
1453 //$ PG 2013-10-17: output of jacobi matrix, IS SLOW, should only be used for debugging reasons
OutputJacobian(const Matrix & M) const1454 void NumNLSolver::OutputJacobian(const Matrix& M) const
1455 {
1456 ostringstream string;
1457 M.PrintToMatlab(string); //extract string from stream via string.str().c_str()
1458 SolverUO() << "Jacobian (# " << jaccount << ") = " << string.str().c_str() << "\n";
1459 }
1460
1461 //$ PG 2013-10-17: output of sparse jacobi matrix, IS SLOW, should only be used for debugging reasons
OutputJacobian(const SparseMatrix & M) const1462 void NumNLSolver::OutputJacobian(const SparseMatrix& M) const
1463 {
1464 SolverUO() << "WARNING: Output of sparse jacobi matrices uses routine for full matrices!\n";
1465 OutputJacobian(M.GetMatrix());
1466 }
1467
1468
1469
1470
1471