1 //#**************************************************************
2 //#
3 //# filename: linalg.h
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 #pragma once
28
29
30 #define gen_sparsemat_off
31
32 int GenVecCnt();
33 int GenMatCnt();
34
35 void GenSparsemat();
36
37 int GetGenSparsemat();
38
39
40 //-----------------------------------------------------------------
41 //-------------- CLASS VECTOR -------------------------------
42 //-----------------------------------------------------------------
43 //Generation of a vector with user-defined length
44 //Operations: +,-,* (also matrices), <<
45 class Vector
46 {
47 public:
48
49 //Constructors, Destructor
Vector()50 Vector(): l(0),vec(NULL),lalloc(0),ownmemory(0) {};
51 Vector(int leni, int nofill=0, double fill=0);
52 Vector(const Vector& veci);
53 Vector(const Vector& veci, int noownmem);
Vector(double x0,double x1,double x2,double x3,double x4,double x5,double x6)54 Vector(double x0, double x1, double x2, double x3, double x4, double x5, double x6)
55 {
56 GenVector(7); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4; vec[5]=x5; vec[6]=x6;
57 }
58
Vector(double x0,double x1,double x2,double x3,double x4,double x5)59 Vector(double x0, double x1, double x2, double x3, double x4, double x5)
60 {
61 GenVector(6); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4; vec[5]=x5;
62 }
Vector(double x0,double x1,double x2,double x3,double x4)63 Vector(double x0, double x1, double x2, double x3, double x4)
64 {
65 GenVector(5); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4;
66 }
Vector(double w,double x,double y,double z)67 Vector(double w, double x, double y, double z)
68 {
69 GenVector(4);vec[0]=w; vec[1]=x; vec[2]=y; vec[3]=z;
70 }
Vector(double x,double y,double z)71 Vector(double x, double y, double z)
72 {
73 GenVector(3);vec[0]=x; vec[1]=y; vec[2]=z;
74 }
Vector(double x,double y)75 Vector(double x, double y)
76 {
77 GenVector(2);vec[0]=x; vec[1]=y;
78 }
Vector(double x)79 Vector(double x)
80 {
81 GenVector(1);vec[0]=x;
82 }
Vector(const TArray<double> & ta)83 Vector(const TArray<double>& ta){SetVector(ta);} //(RL)
84
85
SetVector(const TArray<double> & ta)86 void SetVector(const TArray<double>& ta) //(RL)
87 {
88 GenVector(ta.Length());
89 for (int i=1; i <= ta.Length(); i++)
90 {
91 vec[i-1]=ta(i);
92 }
93 }
94
95 Vector(const char* vector, char osbc = '[', char csbc = ']', char commac = ',');
96
Vector(Vector3D & v1,Vector3D & v2)97 Vector(Vector3D& v1, Vector3D& v2)
98 {
99 GenVector(6); vec[0]=v1.X(); vec[1]=v1.Y(); vec[2]=v1.Z(); vec[3]=v2.X(); vec[4]=v2.Y(); vec[5]=v2.Z();
100 }
~Vector()101 virtual ~Vector()
102 {
103 if (vec!=NULL && ownmemory && lalloc)
104 {delete [] vec; vec=NULL;}
105 };
106
Init()107 virtual void Init()
108 {
109 l = 0;
110 vec = NULL;
111 lalloc = 0;
112 ownmemory = 0;
113 };
114 //Generates a vector with length "leni" (leni=0 possible)
115 virtual void GenVector(int leni);
116 //Link with other vector, no own memory!!!
LinkWith(const Vector & v)117 virtual void LinkWith(const Vector& v)
118 {
119 if (vec!=NULL && ownmemory && lalloc)
120 {delete [] vec; vec=NULL;}
121 ownmemory = 0;
122 vec = v.vec;
123 l = v.l;
124 lalloc = 0;
125 }
LinkWith(double * ptr,int len)126 virtual void LinkWith(double* ptr, int len)
127 {
128 if (vec!=NULL && ownmemory && lalloc)
129 {delete [] vec; vec=NULL;}
130 ownmemory = 0;
131 vec = ptr;
132 l = len;
133 lalloc = 0;
134 }
135 //Link with other vector, no own memory!!!
LinkWith(const Vector & v,int pos,int newlen)136 virtual void LinkWith(const Vector& v, int pos, int newlen)
137 {
138 if (vec!=NULL && ownmemory && lalloc) {delete [] vec; vec=NULL;}
139 ownmemory = 0;
140 vec = v.vec+(pos-1); //no sizeof(double) !!!!!!!!!!!!
141 l = newlen;
142 lalloc = 0;
143 if (newlen == 0) vec = NULL;
144 }
145
146 virtual int IsValid(double maxval);
147
SetAll(double x)148 virtual void SetAll(double x)
149 {
150 for(int i=0; i< l; i++) vec[i]=x;
151 }
152
SetAll(double x,int i1,int i2)153 virtual void SetAll(double x, int i1, int i2)
154 {
155 for(int i = i1-1; i <= i2-1; i++) vec[i]=x;
156 }
157
SetRandom()158 virtual void SetRandom()
159 {
160 for (int i=0; i<l; i++) vec[i] = (double)rand()/(double)RAND_MAX;
161 }
162
IsZero(double eps)163 virtual int IsZero(double eps) const //check if it is smaller than tolerance eps
164 {
165 int found = 0;
166 int i = 0;
167 if (eps == 0)
168 {
169 while (!found && i < l)
170 {
171 if (vec[i++] != 0) found = 1;
172 }
173 }
174 else
175 {
176 while (!found && i < l)
177 {
178 if (fabs(vec[i++]) > eps) found = 1;
179 }
180 }
181 return !found;
182 }
183
Mult(double x)184 virtual void Mult(double x)
185 {
186 for(int i=0; i< l; i++) vec[i]*=x;
187 }
188
189 virtual void SetLen(int leni);
190
191 //Generate vector with 3 components
Set3D(double x,double y,double z)192 virtual void Set3D(double x,double y,double z)
193 {
194 if (l!=3)
195 {
196 if (vec!=NULL && ownmemory && lalloc) {delete[] vec;}
197 GenVector(3);
198 }
199 vec[0]=x; vec[1]=y; vec[2]=z;
200 };
201
202 //Generate vector with 2 components
Set2D(double x,double y)203 virtual void Set2D(double x,double y)
204 {
205 if (l!=2)
206 {
207 if (vec!=NULL && ownmemory && lalloc) {delete[] vec;}
208 GenVector(2);
209 }
210 vec[0]=x; vec[1]=y;
211 };
212
213 //Generate vector with 6 components
Set6D(double v1,double v2,double v3,double v4,double v5,double v6)214 virtual void Set6D(double v1, double v2, double v3, double v4, double v5, double v6)
215 {
216 if (l!=6)
217 {
218 if (vec!=NULL && ownmemory && lalloc) {delete[] vec;}
219 GenVector(6);
220 }
221 vec[0]=v1;
222 vec[1]=v2;
223 vec[2]=v3;
224 vec[3]=v4;
225 vec[4]=v5;
226 vec[5]=v6;
227 };
228
229 //Assignment-operator
230 virtual Vector& operator= (const Vector& veci);
231
232 //Special Assignment-operator
233 virtual Vector& operator= (const Vector3D& veci);
234
235 //Logical operator
236 virtual int operator== (const Vector& vec1);
237
238
239
240 //Referencing access-operator, ZERO-based
241 virtual double& operator[](int elem)
242 {
243 #ifndef __QUICKMATH
244 release_assert((elem>=0) && (elem<l));
245 #endif
246 return vec[elem];
247 };
248
249
250 //Referencing ccess-operator for constant access, ZERO-based
251 virtual const double& operator[](int elem) const
252 {
253 #ifndef __QUICKMATH
254 release_assert((elem>=0) && (elem<l));
255 #endif
256 return vec[elem];
257 };
258
259
260 //Referencing access-operator
operator()261 virtual double& operator()(int elem)
262 {
263 #ifndef __QUICKMATH
264 release_assert((elem>0) && (elem<=l));
265 #endif
266 return vec[elem-1];
267 };
268
269
270 //Referencing access-operator for constant access
operator()271 virtual const double& operator()(int elem) const
272 {
273 #ifndef __QUICKMATH
274 release_assert((elem>0) && (elem<=l));
275 #endif
276 return vec[elem-1];
277 };
278
279
280 //Returns the elem-component of a vector
Elem(int elem)281 virtual double& Elem(int elem)
282 {
283 #ifndef __QUICKMATH
284 release_assert((elem>0) && (elem<=l));
285 #endif
286 return vec[elem-1];
287 };
288
289 //Returns the constant elem-component of a vector
Get(int elem)290 virtual const double& Get(int elem) const
291 {
292 #ifndef __QUICKMATH
293 release_assert((elem>0) && (elem<=l));
294 #endif
295 return vec[elem-1];
296 };
297
298 //Returns the elem-component of a vector, 0-based
Elem0(int elem)299 virtual double& Elem0(int elem)
300 {
301 #ifndef __QUICKMATH
302 release_assert((elem>=0) && (elem<l));
303 #endif
304 return vec[elem];
305 };
306
307 //Returns the constant elem-component of a vector, 0-based
Get0(int elem)308 virtual const double& Get0(int elem) const
309 {
310 #ifndef __QUICKMATH
311 release_assert((elem>=0) && (elem<l));
312 #endif
313 return vec[elem];
314 };
315
GetVecPtr()316 virtual double* GetVecPtr() const {return vec;}
GetVecPtr(int i)317 virtual double* GetVecPtr(int i) const {return &(vec[i-1]);}
318
319 //Returns the length of a vector
GetLen()320 virtual int GetLen() const {return l; };
Length()321 virtual int Length() const {return l; };
322
323 //Returns the quadratic norm of a vector
324 virtual double GetNorm() const;
325
326 //Returns the maximum-norm of a vector
327 virtual double MaxNorm() const;
328
329 //Returns the minimum-norm of a vector
330 virtual double MinNorm() const;
331
332 //Returns the normalized planar perpendicular vector of a vector
333 virtual Vector GetNormed2DNormalVec() const;
334
335 //Returns the normalized planar vector of a vector
336 virtual Vector GetNormed2DVec() const;
337
338 //Returns the cross sum of the components of a vector
339 virtual double Sum() const;
340
341 //Returns the cross sum of the first n components of a vector
342 virtual double Sum(int n) const;
343
344 //Fills a vector with zeros
345 virtual void FillWithZeros();
346
347 //Inserts a vector v at reference point ref
348 // void FillInVector(const Vector& v, const IVector& ref);
349
350 //Inserts a vector v at reference point ref and returns it
351 virtual Vector Append(const Vector& v) const;
352
353 virtual void Insert(const Vector& v, int pos);
Copy(const Vector & v,int vpos,int thispos,int len)354 virtual void Copy(const Vector& v, int vpos, int thispos, int len)
355 {
356 int i;
357 for (i = 0; i < len; i++)
358 {
359 Elem(i+thispos) = v(i+vpos);
360 }
361 }
362
363 //Returns a fraction of a vector
364 virtual Vector SubVector(int from, int to) const;
365
366 //add k*v to vector
367 virtual void MultAdd(const double& k, const Vector& v);
368
369 //Arithmetic operations with one parameter
370 virtual Vector& operator+= (const Vector& v1);
371 virtual Vector& operator-= (const Vector& v1);
372 virtual Vector& operator*= (const double& v);
373
374 //Arithmetic operations with two parameters
375 friend Vector operator+ (const Vector& v1, const Vector& v2);
376 friend Vector operator- (const Vector& v1, const Vector& v2);
377 friend Vector operator* (const Vector& vec, const double& val);
378 friend Vector operator* (const double& val, const Vector& vec);
379 friend Vector operator* (const Vector& v, const Matrix& m);
380 friend Vector operator* (const Matrix& m, const Vector& v);
381 friend double operator* (const Vector& vec1, const Vector& vec2);
382 friend Vector3D operator* (const Matrix3D& m, const Vector& v);
383 friend Vector2D operator* (const Matrix2D& m, const Vector& v);
384 friend void Mult(const Matrix2D& m, const Vector2D& v, Vector& res);
385
386 friend void Mult(const Matrix& m, const Vector& v, Vector& res); //computes res=m*v
387 friend void MultTp(const Matrix& m, const Vector& v, Vector& res); //computes res=m.GetTp()*v
388 friend void Mult(const Matrix& m, const Vector& v, Vector3D& res); //computes res=m*v
389 friend void MultTp(const Matrix& m, const Vector& v, Vector3D& res); //computes res=m*v
390 friend void Mult(const Vector& v, const Matrix& m, Vector& res); //computes res=m*v
391 friend void Mult(const Vector& v1, const Vector& v2, Matrix& res); //computes res=v1*v2^T, where v1 and v2 are interpreted as a column vectors
392 friend void MultBW(const Matrix& m, const Vector& v, Vector& res, int bw); //computes res=m*v with bandwidth bw
393 friend void Mult(const Matrix3D& m, const Vector3D& v, Vector& res);
394 friend void Mult(const MatrixXD& m, const Vector& v, Vector& res);
395 friend void Mult(const Matrix& m, const Vector3D& v, Vector& res); //computes res=m*v
396 friend void Mult(const Matrix& m, const Vector2D& v, Vector& res); //computes res=m*v
397
398 friend void Mult(const SparseMatrix& m, const Vector& v, Vector& res); //computes res=m*v
399
400 //transform strain and stress vectors to matrices and vice versa (by PG)
401 friend void StrainVectorToMatrix2D(Matrix2D& m, const Vector& v);
402 friend void StrainVectorToMatrix2D(Matrix2D& m, const Vector3D& v);
403 friend void StressVectorToMatrix2D(Matrix2D& m, const Vector& v);
404 friend void Matrix2DToStrainVector(Vector& v, const Matrix2D& m);
405 friend void Matrix2DToStressVector(Vector& v, const Matrix2D& m);
406
407 friend void StrainVectorToMatrix3D(Matrix3D& m, const Vector& v);
408 friend void StressVectorToMatrix3D(Matrix3D& m, const Vector& v);
409 friend void Matrix3DToStrainVector(Vector& v, const Matrix3D& m);
410 friend void Matrix3DToStressVector(Vector& v, const Matrix3D& m);
411
412
413 // Vector entries are equally distributed values of piecewise linear function on [-1, 1]
414 // Interpolate returns value of function in point ploc � [-1, 1]
415 // Vector entries specify function values, ordered as
416 // [ 1 2 3 ... width ] --> [x=-1 x=1]
417 virtual double Interpolate(const double ploc) const;
418
419 //Output parameter
420 mystr& MakeString(mystr intro = mystr(""));
421 friend ostream& operator<<(ostream& os, const Vector& v);
422
423 protected:
424
425 double* vec;
426 int l;
427 int lalloc;
428 int ownmemory;
429 };
430
431 template <int data_size>
432 class ConstVector: public Vector
433 {
434 public:
435
436 //Constructors, Destructor
ConstVector()437 ConstVector()
438 {
439 lalloc = 0;
440 ownmemory = 0;
441 l = data_size;
442 vec = &constdata[0];
443 SetAll(0.);
444 }
ConstVector(const ConstVector & veci)445 ConstVector(const ConstVector& veci)
446 {
447 lalloc = 0;
448 ownmemory = 0;
449 vec = &constdata[0];
450
451 #ifndef __QUICKMATH
452 release_assert(data_size >= veci.l);
453 #endif
454 l = veci.l;
455 for (int i=0; i < l; i++)
456 {
457 constdata[i] = veci.vec[i];
458 }
459 }
460
461 ConstVector(int leni, int nofill=0, double fill=0)
462 {
463 lalloc = 0;
464 ownmemory = 0;
465 l = leni;
466 vec = &constdata[0];
467 if (!nofill) SetAll(fill);
468 }
469
ConstVector(double x0)470 ConstVector(double x0)
471 {
472 release_assert(data_size >= 1);
473 InitConstVector(); vec[0]=x0;
474 }
ConstVector(double x0,double x1)475 ConstVector(double x0, double x1)
476 {
477 release_assert(data_size >= 2);
478 InitConstVector(); vec[0]=x0; vec[1]=x1;
479 }
ConstVector(double x0,double x1,double x2)480 ConstVector(double x0, double x1, double x2)
481 {
482 release_assert(data_size >= 3);
483 InitConstVector(); vec[0]=x0; vec[1]=x1; vec[2]=x2;
484 }
ConstVector(double x0,double x1,double x2,double x3)485 ConstVector(double x0, double x1, double x2, double x3)
486 {
487 release_assert(data_size >= 4);
488 InitConstVector(); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3;
489 }
ConstVector(double x0,double x1,double x2,double x3,double x4)490 ConstVector(double x0, double x1, double x2, double x3, double x4)
491 {
492 release_assert(data_size >= 5);
493 InitConstVector(); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4;
494 }
ConstVector(double x0,double x1,double x2,double x3,double x4,double x5)495 ConstVector(double x0, double x1, double x2, double x3, double x4, double x5)
496 {
497 release_assert(data_size >= 6);
498 InitConstVector(); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4; vec[5]=x5;
499 }
ConstVector(double x0,double x1,double x2,double x3,double x4,double x5,double x6,double x7,double x8,double x9)500 ConstVector(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9)
501 {
502 release_assert(data_size >= 10);
503 InitConstVector(); vec[0]=x0; vec[1]=x1; vec[2]=x2; vec[3]=x3; vec[4]=x4; vec[5]=x5; vec[6]=x6; vec[7]=x7; vec[8]=x8; vec[9]=x9;
504 }
InitConstVector()505 virtual void InitConstVector()
506 {
507 lalloc = 0;
508 ownmemory = 0;
509 l = data_size;
510 vec = &constdata[0];
511 }
GetDataSize()512 virtual int GetDataSize() const {return data_size;}
513
514
515 virtual ConstVector<data_size> & operator= (const Vector& veci)
516 {
517 if (this == &veci) {return *this;}
518 if (veci.Length()!=0)
519 {
520 #ifndef __QUICKMATH
521 release_assert(data_size>=veci.Length());
522 #endif
523 l = veci.Length();
524 memcpy(vec,&veci[0],l*sizeof(double));
525 }
526 else
527 {
528 this->l=0;
529 //lalloc = 0; vec = NULL;
530 }
531 return *this;
532 }
533
534 private:
535 double constdata[data_size];
536 };
537
538
539
540
541 //-----------------------------------------------------------------
542 //-------------- CLASS MATRIX -------------------------------
543 //-----------------------------------------------------------------
544 //Generation of a Matrix with user-defined size
545 //Operations: +,-,*,<<
546 class Matrix
547 {
548 public:
549
550 //Constructors,destructor
Matrix()551 Matrix()
552 {
553 rows=0;
554 cols=0;
555 mat=NULL;
556 lalloc = 0;
557 };
558
559 Matrix(int rowsi, int colsi, int nofill=0, double fill=0);
560 Matrix(double x);
561 Matrix(double x11, double x12, double x21, double x22);
562 Matrix(double x11, double x12, double x13,
563 double x21, double x22, double x23,
564 double x31, double x32, double x33);
565
566 Matrix(const Matrix& mati);
567 Matrix(const Matrix3D mat3d);
568
~Matrix()569 virtual ~Matrix()
570 {
571 if (mat!=NULL && lalloc != 0) {delete [] mat; mat=NULL;}
572 };
573
574
Matrix(int r,int c,double * ptr)575 Matrix(int r, int c, double * ptr)
576 {
577 lalloc = 0;
578 mat = ptr;
579 rows = r;
580 cols = c;
581 }
582
Init()583 virtual void Init()
584 {
585 rows=0;
586 cols=0;
587 mat=NULL;
588 lalloc = 0;
589 };
IsConstSizeMatrix()590 virtual int IsConstSizeMatrix() const {return 0;}
591
592 //Setsize of matrix, don't change if ok
SetSize(int rowsi,int colsi)593 virtual void SetSize(int rowsi, int colsi)
594 {
595 if (rows!=rowsi || cols!=colsi || mat == NULL) {Resize(rowsi,colsi,1);}
596 }
GetMatPtr()597 virtual double* GetMatPtr() const {return mat;}
MaxAllocatedSize()598 virtual int MaxAllocatedSize() const {return lalloc;}
599
600 //Set new size of matrix, delete old one
601 virtual void Resize(int rowsi, int colsi, int nofill=0, double fill=0);
602
603 virtual void CopyFrom(const Matrix& m, int row1, int col1, int row2, int col2);
604
605 virtual void CopyFrom(const Matrix& m, int row1, int col1, int row2, int col2, const IVector& r);
606
607 virtual void CopyFrom(const SparseMatrix& m, int row1, int col1, int row2, int col2);
608
609 virtual void CopyFrom(const SparseMatrix& m, int row1, int col1, int row2, int col2, const IVector& r);
610
611 //Access-functions for private data
Getrows()612 virtual int Getrows() const {return rows; };
Getcols()613 virtual int Getcols() const {return cols; };
614
615 //Assignment-operator
616 virtual Matrix& operator= (const Matrix& mati);
617
618 //Assignment-operator
619 virtual Matrix& operator= (const double& val)
620 {
621 for (int i=0; i < rows*cols; i++)
622 {
623 mat[i] = val;
624 }
625 return *this;
626 }
627
628 //Returns Determinate of matrix
629 virtual double Det() const;
630
631 //Returns the transposed matrix
632 virtual Matrix GetTp() const;
633
634
635 virtual double* operator[](int row0) //zero based for [][] access
636 {
637 #ifndef __QUICKMATH
638 release_assert((row0>=0) && (row0<rows));
639 #endif
640 return &(mat[row0*cols]);
641 };
642
643 virtual const double* operator[](int row0) const //zero based for [][] access
644 {
645 #ifndef __QUICKMATH
646 release_assert((row0>=0) && (row0<rows));
647 #endif
648 return &(mat[row0*cols]);
649 };
650
651 /*
652 //Referencing access-operator on element using absolute indexing
653 double& operator[](int elem)
654 {
655 #ifndef __QUICKMATH
656 release_assert((elem>0) && (elem<=rows*cols));
657 #endif
658 return mat[elem-1];
659 };
660
661
662 //Referencing constant access-operator on element using absolute indexing
663 const double& operator[](int elem) const
664 {
665 #ifndef __QUICKMATH
666 release_assert((elem>0) && (elem<=rows*cols));
667 #endif
668 return mat[elem-1];
669 };
670 */
671
672
673 //Referencing access-operator on element using row- and column-values
operator()674 virtual double& operator()(int row, int col)
675 {
676 #ifndef __QUICKMATH
677 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
678 #endif
679 return mat[(row-1)*cols+col-1];
680 };
681
682
683 //Referencing constant access-operator on element using row- and column-values
operator()684 virtual const double& operator()(int row, int col) const
685 {
686 #ifndef __QUICKMATH
687 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
688 #endif
689 return mat[(row-1)*cols+col-1];
690 };
691
692
693 //Referencing access-operator on element using row- and column-values
Elem(int row,int col)694 virtual double& Elem(int row, int col)
695 {
696 #ifndef __QUICKMATH
697 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
698 #endif
699 return mat[(row-1)*cols+col-1];
700 };
701
702
703 //Referencing constant access-operator on element using row- and column-values
Get(int row,int col)704 virtual const double& Get(int row, int col) const
705 {
706 #ifndef __QUICKMATH
707 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
708 #endif
709 return mat[(row-1)*cols+col-1];
710 };
711
712 //Referencing access-operator on element using row- and column-values, 0-based
Get0(int row,int col)713 virtual const double& Get0(int row, int col) const
714 {
715 #ifndef __QUICKMATH
716 release_assert((row>=0) && (row<rows) && (col>=0) && (col<cols));
717 #endif
718 return mat[row*cols+col];
719 };
720
721 //Referencing access-operator on element using row- and column-values, 0-based
Elem0(int row,int col)722 virtual double& Elem0(int row, int col)
723 {
724 #ifndef __QUICKMATH
725 release_assert((row>=0) && (row<rows) && (col>=0) && (col<cols));
726 #endif
727 return mat[row*cols+col];
728 };
729
730
LinkWith(int r,int c,double * ptr)731 virtual void LinkWith(int r, int c, double * ptr)
732 {
733 if (mat!=NULL && lalloc != 0)
734 {delete [] mat; mat=NULL;}
735 lalloc = 0;
736 mat = ptr;
737 rows = r;
738 cols = c;
739 }
740
741 //Access-operator for symmetric band matrix, returns value of element
BGet(int row,int col)742 virtual double BGet(int row, int col) const
743 {
744 if (row > col) {Swap (row,col);}
745 release_assert((row>0) && (row<=rows) && (col>0) && (col<=rows));
746 if ((col>row+cols-1)) {return 0;}
747 return mat[(row-1)*cols+col-row];
748 };
749
750
751 //Access-operator for symmetric band matrix, sets value of element
Bset(int row,int col,const double & v)752 virtual void Bset(int row, int col, const double& v)
753 {
754 if (row > col) {Swap (row,col);}
755 release_assert((row>0) && (row<=rows) && (col<row+cols));
756 if ((col>row+cols-1)) {return;}
757 mat[(row-1)*cols+col-row]=v;
758 };
759
760
761 //Converts symmetric band-matrix to matrix
762 virtual Matrix BToMatrix();
763
764 //Converts symmetric matrix to symmetric band-matrix
765 virtual Matrix MatrixToB();
766
767 //Fills matrix with zeros
768 virtual void FillWithZeros();
769
770 virtual void SetAll(double x);
771
772 //Fills matrix with zeros
773 virtual void FillWithZeros(int row, int col, int nrows, int ncols);
774
775 //Returns the maximum-norm (largest value in matrix)
776 virtual double MaxNorm() const;
777
778
779 //Returns the minimum-norm (smallest value in matrix)
780 virtual double MinNorm() const;
781
782 //Returns the quadratic-norm
783 virtual double Norm2() const;
784
785 //Returns true if matrix is symmetric
786 virtual int IsSymmetric();
787 virtual int IsSymmetric(double eps); //version with small tolerance epsilon
788
789 virtual void MakeSymmetric(); //makes a matrix really symmetric (if there are small errors)
790
791 //Returns true if matrix is quadratic
IsSquare()792 virtual int IsSquare() const {return rows==cols;}
793
794
795 //Returns the bandwidth of a matrix
796 virtual int GetBandWidth() const;
797
GetBandMatrix(Matrix & bm,int lu)798 virtual void GetBandMatrix(Matrix& bm, int lu)
799 {
800 int col_band = lu*2+1;
801 int n = rows;
802
803 bm.SetSize(n,col_band);
804
805 int ind;
806 for (int i = 1; i <= n; i++)
807 {
808 for (int j = 1; j <= col_band; j++)
809 {
810 ind = i-(lu+1)+j;
811 if (ind >= 1 && ind <= n)
812 {
813 bm(i,j) = Get(i, ind);
814 }
815 else bm(i,j) = 0;
816 }
817 }
818 }
819
820 //Returns the column-vector at clomn col
821 virtual Vector GetColVec(int col) const;
822 virtual void GetColVec(int col, Vector& v) const;
823
824 //Returns the row-vector at row row
825 virtual Vector GetRowVec(int row) const;
826 virtual void GetRowVec(int row, Vector& v) const;
827
828
829 //Returns the diagonal values
830 virtual Vector GetDiag() const;
831
832 //Set Matrix to diagonal matrix with value v
SetDiagMatrix(double v)833 virtual void SetDiagMatrix(double v)
834 {
835 #ifndef __QUICKMATH
836 release_assert((rows==cols));
837 #endif
838 for (int i=1; i <= rows; i++)
839 {
840 for (int j=1; j <= rows; j++)
841 {
842 if (i!=j) Elem(i,j) = 0;
843 else Elem(i,j) = v;
844 }
845 }
846 }
847
848 //Returns the Trace
849 virtual double Trace() const;
850
851 //Sets the column-vector at clomn col
852 virtual void SetColVec(const Vector& x, int col);
853
854 //Sets the row-vector at row row
855 virtual void SetRowVec(const Vector& x, int row);
856
857 //Transposes a matrix
858 virtual void TpYs();
859
Mult(double x)860 virtual void Mult(double x)
861 {
862 for(int i=0; i< rows*cols; i++) mat[i]*=x;
863 }
864
865 //Take Submatrix sm at (sr,sr), multiply with x and add it to local matrix at (r,c), size: (nr x nc)
866 virtual void AddSubmatrix(const Matrix& sm, int sr, int sc, int r, int c, int nr, int nc, double x);
867 //Set whole Submatrix sm at *this matrix at (r,c)
868 virtual void SetSubmatrix(const Matrix& sm, int r, int c);
869 virtual void SetSubmatrix(const Matrix& sm, int r, int c, double x);
870 virtual void SetSubmatrix(const Matrix3D& sm, int r, int c);
871
872 //Multiplies the cloumn col with val
873 virtual void MulCol(int col, double val);
874
875 //Multiplies the row row with val
876 virtual void MulRow(int row, double val);
877
878
879 //Insert (smaller) Matrix at position (row,col), delete entries
880 virtual void InsertMatrix(int row, int col, const Matrix& m);
881
882 //Insert (smaller) Matrix at position (row,col), mrows and mcols of matrix m, delete entries
883 virtual void InsertMatrix(int row, int col, const Matrix& m, int fromrow, int fromcol, int mrows, int mcols);
884
885 //Add (possibly smaller) Matrix at position (row,col)
886 virtual void AddMatrix(int row, int col, const Matrix& m);
887 //Add (possibly smaller) Matrix at position given by references
888 virtual void AddMatrix(const TArray<int>& rowref, const TArray<int>& colref, const Matrix& m);
889
890 //Adds the vector vec to row row
891 //Adds the vector vec to row row
892 virtual void AddRowVec(int row, const Vector& vec);
893
894 virtual void AddColVec(int col, const Vector& vec);
895
896 //Adds fromrow multiplied with fact to row row
897 virtual void AddRowVec(int fromRow, int toRow, double fact);
898 //Adds fromrow multiplied with fact to row row
899 virtual void AddRowVec(int fromRow, int toRow, double fact, int fromCol, int toCol);
900
901 //Multiply column col with column colm of matrix m
902 virtual double MultCol(const Matrix& m, int mcol, int col);
903
904 //Adds fromrow multiplied with fact to row row
905 virtual void AddColVec(int fromCol, int toCol, double fact);
906
907 //Sets the column-vector at clomn col, taken from Matrix m, column mcol
908 virtual void SetColVec(const Matrix& m, int mcol, int col);
909
910
911 //Adds fromrow multiplied with fact to row row, including bandwidth
912 virtual void AddRowVecB(int fromRow, int toRow, const double& fact, int bw);
913
914 //Swaps the rows r1 and r2
915 virtual void SwapRows(int r1, int r2);
916
917 //Swaps the columns c1 and c2
918 virtual void SwapCols(int c1, int c2);
919
920 //Multiplies band-matrix with vector vec
921 virtual Vector Bmult(Vector& v) const;
922
923
924 //Inverts matrix, returns 1 if successfull
925 virtual int Invert();
926
927 //Inverts matrix, returns 1 if successfull, faster than invert()
928 virtual int Invert2();
929
930 //Inverts matrix, returns 1 if successfull, using LAPACK
931 virtual int InvertLapack();
932
933 //Computes an estimate of the reciprocal of the condition number of this matrix, returns 1 if successfull, using LAPACK
934 virtual int EstimateReciprocalConditionNumber(double& rcond) const;
935
936 //$ PG 2013-10-30: used for solving an undetermined system (underdetermined: minimum norm solution |q|->min subject to Aq=f; overdetermined: least squares solution |Aq-f|^2->min)
937 // f (input) right hand side
938 // q (output) solution
939 // rcond (input) reciprocal of estimated condition number of A (this matrix), used for recognizing linearly dependendt equations numerically
940 // rank (output) computed rank of matrix A, depending on rcond.
941 virtual int UndeterminedSystemSolve(const Vector& f, Vector& q, double rcond, int& rank);
942
943 //Solving equation mat*q=f using gaussian elimination method, returns 1 if successfull
944 virtual int Solve(const Vector& fv, Vector& q); //K*q=f-->q
945
946 //Solving equation mat*q=f using LAPACK, returns 1 if successfull
947 virtual int SolveLapack(Vector& q);
948
949 //LU-Decomposition of matrix, returns 1 if successfull, indx is vector of row interchanges
950 virtual int LU(IVector& indx);
951
952 //Solve A*x=b, A is decomposed with LU, result is x=b
953 virtual int LUBCKSUB(IVector& indx, Vector& b);
954
955 //Solve A*x=b, A is decomposed with LU, result is x=b
956 //virtual int LUBCKSUB(IVector& indx, SparseVector& b, const Vector& bdense);
957
958 //Solving equation mat*q=f using gaussian elimination method
959 //and considering band-structure, returns 1 if successfull
960 virtual int Bsolve(Vector f, Vector& q); //K*q=f-->q
961
962
963 //Solving the equation mat*q=f using a band-matrix based on Cholesky-Method, returns 1 if successfull
964 virtual int CholeskyBsolve(Vector f, Vector& q);
965
966 virtual int QRDecomposition(Matrix& Q, Matrix& R);
967
968 virtual int Eigenvalues(Vector& ev);
969 virtual int EigenvaluesLapack(Vector& lami) const;
970
971 //writes out matrix using MatLab format
972 virtual void PrintToMatlab(ostream& os) const;
973 virtual void PrintToMatlabSparse(ostream& os, double tol=1e-16);
974
975 //writes out matrix using Maple format
976 virtual void PrintToMaple(ostream& os);
977
978 //writes out matrix using Mathematica (Table) format
979 virtual mystr PrintToMathematica();
980
981
982
983
984 //Biconjugate gradient method
985 virtual int BCG(Vector f, Vector& q, const Vector& qs);
986
987 // Matrix entries are equally distributed values of piecewise bilinear function on the unit square
988 // Interpolate returns value of function in point ploc � [-1, 1] x [-1, 1]
989 // Matrix entries specify function values, ordered as
990 // [ 1 2 3 ... width ] [ y=1 ]
991 // [ w+1 w+2 w+3 ... 2w ] -----> [x=-1 x=1]
992 //
993 // [ (h-1)w+1 (h-1)w+1 (h-1)w+2 ... hw] [ y=-1 ]
994 virtual double Interpolate(const Vector2D& ploc) const;
995 // Interpolate at given x-coord xloc,
996 // y-coord corresponds to vector entries (at exit, vector length = number of rows of matrix
997 void InterpolateX(const double xloc, Vector& interpol_vector) const;
998
999 //Arithmetic operations with 1 parameter
1000 virtual Matrix& operator+= (const Matrix& m);
1001 virtual Matrix& operator-= (const Matrix& m);
1002 virtual Matrix& operator*= (const double& val);
1003
1004 //$ SW 2013-10-9: added
1005 //calculates *= A where A is a square matrix (e.g. a rotation matrix) of dimension 3x3, 2x2 or 1x1
1006 virtual void ApplySqrMat(const MatrixXD &A);
1007
1008 // Write Matrix to string to export to C "double intro[rows][cols] = {.,.,.};"
1009 mystr& MakeString(mystr intro = mystr(""));
1010
1011 //Arithmetic operations with 2 parameters
1012 friend Matrix operator+ (const Matrix& m1, const Matrix& m2);
1013 friend Matrix operator- (const Matrix& m1, const Matrix& m2);
1014 friend Matrix operator* (const Matrix& mat, const double& val);
1015 friend Matrix operator* (const double& val, const Matrix& mat);
1016 friend Matrix operator* (const Matrix& m1, const Matrix& m2);
1017 friend Vector operator* (const Vector& v, const Matrix& m);
1018 friend Vector operator* (const Matrix& m, const Vector& v);
1019 friend void Mult(const Matrix& m, const Vector& v, Vector& res); //computes res=m*v
1020 friend void Mult(const Matrix& m, const Vector& v, Vector3D& res); //computes res=m*v
1021 friend void MultTp(const Matrix& m, const Vector& v, Vector3D& res); //computes res=m*v
1022 friend void Mult(const Vector& v, const Matrix& m, Vector& res); //computes res=m*v
1023 friend void Mult(const Vector& v1, const Vector& v2, Matrix& res); //computes res=v1*v2^T, where v1 and v2 are interpreted as a column vectors
1024 friend void MultBW(const Matrix& m, const Vector& v, Vector& res, int bw); //computes res=m*v with bandwidth bw
1025 friend void Mult(const Matrix& m, const Vector3D& v, Vector& res); //computes res=m*v
1026 friend void Mult(const Matrix& m, const Vector2D& v, Vector& res); //computes res=m*v
1027 friend void Mult(const Matrix& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2
1028 friend void MultTp(const Matrix& m1, const Matrix& m2, Matrix& res); //computes res=m1^T*m2
1029 friend void MultSym(const Matrix& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2, result is assumed to be symmetric
1030 friend void MultSymTp(const Matrix& m1, const Matrix& m2, Matrix& res); //computes res=m1^T*m2, result is assumed to be symmetric
1031 friend void Mult(const MatrixXD& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2
1032 friend void Mult(const Matrix& m1, const MatrixXD& m2, Matrix& res); //computes res=m1*m2
1033 friend Matrix operator* (const SparseMatrix& m1, const Matrix& m2);
1034 friend void Mult(const SparseMatrix& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2
1035
1036 //writes out matrix with constandt width and factor normalized
1037 friend ostream& operator<<(ostream& os, const Matrix& m);
1038
SetMatrix2n(TArray<Vector2D> & list)1039 virtual void SetMatrix2n(TArray<Vector2D>& list) // $ MSax 2013-07-09 : added
1040 {
1041 SetSize(list.Length(),2);
1042 Vector2D tmp;
1043 for (int i=1; i <= list.Length(); i++)
1044 {
1045 tmp = list.Get(i);
1046 mat[2*i-2]=tmp.X();
1047 mat[2*i-1]=tmp.Y();
1048 }
1049 }
1050
GetVector2DList(TArray<Vector2D> & list)1051 virtual void GetVector2DList(TArray<Vector2D>& list) // $ MSax 2013-07-09 : added
1052 {
1053 list.Init();
1054 for (int i=1; i <= rows; i++)
1055 {
1056 list.Add(Vector2D(mat[2*i-2],mat[2*i-1]));
1057 }
1058 }
1059
1060 protected:
1061
1062 double* mat;
1063 int rows, cols;
1064 int lalloc;
1065 };
1066
1067 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1068 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1069 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1070
1071
1072 template <int data_size>
1073 class ConstMatrix: public Matrix
1074 {
1075 public:
1076
1077 //Constructors, Destructor
ConstMatrix()1078 ConstMatrix()
1079 {
1080 InitConstMatrix();
1081 }
1082
ConstMatrix(const ConstMatrix & mati)1083 ConstMatrix(const ConstMatrix& mati)
1084 {
1085 InitConstMatrix();
1086
1087 #ifndef __QUICKMATH
1088 release_assert(data_size >= mati.rows*mati.cols);
1089 #endif
1090 rows = mati.rows;
1091 cols = mati.cols;
1092 for (int i=0; i < rows*cols; i++)
1093 {
1094 constdata[i] = mati.mat[i];
1095 }
1096 }
1097
MaxAllocatedSize()1098 virtual int MaxAllocatedSize() const {return data_size;}
IsConstSizeMatrix()1099 virtual int IsConstSizeMatrix() const {return 1;}
1100
1101 ConstMatrix(int rowsi, int colsi, int nofill=0, double fill=0)
1102 {
1103 #ifndef __QUICKMATH
1104 release_assert(data_size >= rowsi*colsi);
1105 #endif
1106 InitConstMatrix();
1107 Resize(rowsi, colsi, nofill, fill);
1108 }
1109
ConstMatrix(double x0)1110 ConstMatrix(double x0)
1111 {
1112 #ifndef __QUICKMATH
1113 release_assert(data_size >= 1);
1114 #endif
1115 InitConstMatrix();
1116
1117 rows=1; cols=1;
1118 mat[0]=x0;
1119 }
1120
ConstMatrix(double x11,double x12,double x21,double x22)1121 ConstMatrix(double x11, double x12, double x21, double x22)
1122 {
1123 #ifndef __QUICKMATH
1124 release_assert(data_size >= 4);
1125 #endif
1126 InitConstMatrix();
1127
1128 rows=2; cols=2;
1129 mat[0]=x11;
1130 mat[1]=x12;
1131 mat[2]=x21;
1132 mat[3]=x22;
1133 }
1134
ConstMatrix(double x11,double x12,double x13,double x21,double x22,double x23,double x31,double x32,double x33)1135 ConstMatrix(double x11, double x12, double x13,
1136 double x21, double x22, double x23,
1137 double x31, double x32, double x33)
1138 {
1139 #ifndef __QUICKMATH
1140 release_assert(data_size >= 9);
1141 #endif
1142 InitConstMatrix();
1143
1144 rows=3; cols=3;
1145 mat[0]=x11;
1146 mat[1]=x12;
1147 mat[2]=x13;
1148 mat[3]=x21;
1149 mat[4]=x22;
1150 mat[5]=x23;
1151 mat[6]=x31;
1152 mat[7]=x32;
1153 mat[8]=x33;
1154 }
1155
InitConstMatrix()1156 virtual void InitConstMatrix()
1157 {
1158 rows=0;
1159 cols=0;
1160 mat = &constdata[0];
1161 lalloc = 0;
1162 }
GetDataSize()1163 virtual int GetDataSize() const {return data_size;}
1164
1165 private:
1166 double constdata[data_size];
1167 };
1168
1169
1170
1171
1172 Matrix GetDiag(int n, double val = 1);
1173 double Mises(const Matrix& m);
1174 void PrintMatrix01(const Matrix& m);
1175
1176
1177
1178 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1179 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1180 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1181
1182
1183 //-----------------------------------------------------------------
1184 //-------------- CLASS SPARSEVECTOR -------------------------------
1185 //-----------------------------------------------------------------
1186 //Generation of a vector with user-defined length
1187 //Operations: +,-,* (also matrices), <<
1188 class SparseVector
1189 {
1190 public:
1191
1192 //Constructors, Destructor
SparseVector()1193 SparseVector(): l(0),vec(NULL),lalloc(0), rowind(NULL), nelem(0), dnull(0.) {};
1194
SparseVector(int leni,int ilalloc)1195 SparseVector(int leni, int ilalloc)
1196 {
1197 GenVector(leni, ilalloc);
1198 };
1199
SparseVector(const SparseVector & veci)1200 SparseVector(const SparseVector& veci)
1201 {
1202 GenVector(veci.l, veci.lalloc);
1203 nelem = veci.nelem;
1204 if (nelem!=0)
1205 {
1206 memcpy(vec,veci.vec,nelem*sizeof(double));
1207 memcpy(rowind,veci.rowind,nelem*sizeof(int));
1208 }
1209 }
1210
~SparseVector()1211 virtual ~SparseVector()
1212 {
1213 Destroy();
1214 };
1215
1216 //Generates a vector with length "leni" (leni=0 possible)
1217 void GenVector(int leni, int ilalloc = 0)
1218 {
1219 dnull = 0.;
1220
1221 l=leni;
1222 nelem = 0;
1223 lalloc = ilalloc;
1224
1225 if (lalloc)
1226 {
1227 vec = new double[lalloc];
1228 rowind = new int[lalloc];
1229 }
1230 else
1231 {
1232 vec = NULL;
1233 rowind = NULL;
1234 }
1235 }
Destroy()1236 void Destroy()
1237 {
1238 if (vec!=NULL && lalloc) {delete [] vec; vec=NULL;}
1239 if (rowind!=NULL && lalloc) {delete [] rowind; rowind=NULL;}
1240 }
1241
1242 //deletes data!!!
1243 void SetLen(int len, int ilalloc = 0)
1244 {
1245 if (ilalloc > lalloc)
1246 {
1247 Destroy();
1248 GenVector(len, ilalloc);
1249 }
1250 else
1251 {
1252 l = len;
1253 }
1254 }
1255
CopyFrom(const Vector & v)1256 void CopyFrom(const Vector& v)
1257 {
1258 int nonzero = 0;
1259 for (int i = 0; i < v.Length(); i++)
1260 if (v[i] != 0.) nonzero++;
1261
1262 SetLen(v.Length(), nonzero);
1263 FillWithZeros();
1264
1265 for (int i = 0; i < v.Length(); i++)
1266 {
1267 if (v[i] != 0.) AddEntry(i+1, v[i]);
1268 }
1269
1270 }
1271
Length()1272 int Length() const {return l;}
GetLen()1273 int GetLen() const {return l;}
1274
NEntries()1275 int NEntries() const {return nelem;}
NAlloc()1276 int NAlloc() const {return lalloc;}
1277
1278 //const Referencing access-operator
operator()1279 const double& operator()(int row) const
1280 {
1281 for (int i=0; i < nelem; i++)
1282 {
1283 if (row == rowind[i]+1) return vec[i];
1284 }
1285 //static double dnull=0;
1286 return dnull;
1287 };
1288
1289 //Referencing access-operator
operator()1290 double& operator()(int row)
1291 {
1292 int i=0;
1293 while (row < rowind[i]+1 && i < nelem)
1294 {
1295 i++;
1296 }
1297 if (row == rowind[i]+1) return vec[i];
1298
1299 //element does not exist:
1300 if (nelem == lalloc) //resize?
1301 {
1302 int oldlalloc = lalloc;
1303 if (lalloc == 0) lalloc = 1;
1304 else lalloc *= 2;
1305
1306 double* new_vec = new double[lalloc];
1307 int* new_rowind = new int[lalloc];
1308
1309 if (oldlalloc != 0)
1310 {
1311 memcpy(new_vec,vec,oldlalloc*sizeof(double));
1312 memcpy(new_rowind,rowind,oldlalloc*sizeof(int));
1313
1314 delete [] vec;
1315 delete [] rowind;
1316 }
1317 vec = new_vec;
1318 rowind = new_rowind;
1319 }
1320
1321 if (i < nelem)
1322 {
1323 for (int j=nelem; j > i; j--)
1324 {
1325 rowind[j] = rowind[j-1];
1326 vec[j] = vec[j-1];
1327 }
1328 }
1329 nelem++;
1330 vec[i] = 0;
1331 rowind[i] = row-1;
1332 return vec[i];
1333 };
1334
1335
1336
1337 //read Entry i:
GetEntry(int i,int & row,double & d)1338 void GetEntry(int i, int& row, double& d) const
1339 {
1340 #ifndef __QUICKMATH
1341 release_assert((i>0) && (i<=nelem));
1342 #endif
1343 row = rowind[i-1]+1;
1344 d = vec[i-1];
1345 }
1346
SetEntry(int i,int row,double d)1347 void SetEntry(int i, int row, double d) const
1348 {
1349 #ifndef __QUICKMATH
1350 release_assert((i>0) && (i<=nelem));
1351 #endif
1352 rowind[i-1] = row-1;
1353 vec[i-1] = d;
1354 }
1355
Entry(int i)1356 double& Entry(int i)
1357 {
1358 #ifndef __QUICKMATH
1359 release_assert((i>0) && (i<=nelem));
1360 #endif
1361 return vec[i-1];
1362 }
1363
Entry(int i)1364 const double& Entry(int i) const
1365 {
1366 #ifndef __QUICKMATH
1367 release_assert((i>0) && (i<=nelem));
1368 #endif
1369 return vec[i-1];
1370 }
1371
1372 //read Entry i:
AddEntry(int row,double d)1373 void AddEntry(int row, double d)
1374 {
1375 if (nelem == lalloc)
1376 {
1377 int oldlalloc = lalloc;
1378 if (lalloc == 0) lalloc = 1;
1379 else lalloc *= 2;
1380
1381 double* new_vec = new double[lalloc];
1382 int* new_rowind = new int[lalloc];
1383
1384 if (oldlalloc != 0)
1385 {
1386 memcpy(new_vec,vec,oldlalloc*sizeof(double));
1387 memcpy(new_rowind,rowind,oldlalloc*sizeof(int));
1388
1389 delete [] vec;
1390 delete [] rowind;
1391 }
1392 vec = new_vec;
1393 rowind = new_rowind;
1394 }
1395 rowind[nelem] = row-1;
1396 vec[nelem] = d;
1397 nelem++;
1398 }
1399
GetRow(int i)1400 int GetRow(int i) const
1401 {
1402 #ifndef __QUICKMATH
1403 release_assert((i>0) && (i<=nelem));
1404 #endif
1405 return rowind[i-1]+1;
1406 }
1407
FillWithZeros()1408 void FillWithZeros()
1409 {
1410 nelem = 0;
1411 }
1412
1413 private:
1414
1415 //as vector:
1416 int l;
1417 double* vec;
1418
1419 //dynamic:
1420 int lalloc;
1421 int nelem;
1422 int* rowind;
1423
1424 double dnull;
1425 };
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438 //-----------------------------------------------------------------
1439 //-------------- CLASS SPARSEMATRIX -------------------------------
1440 //-----------------------------------------------------------------
1441 //Generation of a Matrix with user-defined size
1442 //Operations: +,-,*,<<
1443
1444 //column based sparse matrix!
1445 //element i, j: find x such that colind[x]==j, mat[rowind[i]+x]
1446 //allocated length of row: through rowind[i+1]-rowind;
1447 //used length of row: rowlen
1448 //column index is stored in colind for every element mat
1449
1450 //undefined element:
1451 #define sparseMatrixVoid -1
1452 //increase factor for reallocate:
1453 #define sparseMatrixInc 2
1454
1455 const double sparsetolzero = 1e-20;
1456
1457 class SparseMatrix
1458 {
1459 public:
1460
1461 //Constructors,destructor
SparseMatrix()1462 SparseMatrix()
1463 {
1464 Init();
1465 };
1466
1467 void GenMat(int lalloc, int rowsi);
1468
Destroy()1469 void Destroy()
1470 {
1471 if (mat!=NULL) {delete [] mat; mat=NULL;}
1472 if (rowind!=NULL) {delete [] rowind; rowind=NULL;}
1473 if (rowlen!=NULL) {delete [] rowlen; rowlen=NULL;}
1474 if (colind!=NULL) {delete [] colind; colind=NULL;}
1475 }
1476
1477 SparseMatrix(int rowsi, int colsi, int initcols=1)
1478 {
1479 if (rowsi*colsi == 0) {Init(); return;}
1480
1481 if (initcols < 1) {initcols = 1;}
1482
1483 rows = rowsi;
1484 cols = colsi;
1485
1486 lalloc = initcols*rowsi;
1487 GenMat(lalloc, rowsi);
1488
1489 for (int i=0; i < rows; i++)
1490 {
1491 rowind[i] = i*initcols;
1492 rowlen[i] = 0;
1493 for (int j=0; j < initcols; j++)
1494 {
1495 colind[rowind[i]+j] = sparseMatrixVoid;
1496 }
1497 }
1498 }
1499
SparseMatrix(const SparseMatrix & mati)1500 SparseMatrix(const SparseMatrix& mati)
1501 {
1502 lalloc = mati.lalloc;
1503 rows = mati.rows;
1504 cols = mati.cols;
1505
1506 GenMat(lalloc, rows);
1507
1508 for (int i = 0; i < rows; i++)
1509 {
1510 rowind[i] = mati.rowind[i];
1511 rowlen[i] = mati.rowlen[i];
1512 }
1513 for (int i = 0; i < lalloc; i++)
1514 {
1515 colind[i] = mati.colind[i];
1516 mat[i] = mati.mat[i];
1517 }
1518 }
1519
~SparseMatrix()1520 virtual ~SparseMatrix()
1521 {
1522 Destroy();
1523 };
1524
SparseMatrix(const Matrix & mati)1525 SparseMatrix(const Matrix& mati)
1526 {
1527 Init();
1528 CopyFrom(mati);
1529 }
1530
1531 void CopyFrom(const Matrix& mati);
1532
1533 void CopyFrom(const Matrix& mati, int row1, int col1, int row2, int col2);
1534
1535 void CopyFrom(const SparseMatrix& mati, int row1, int col1, int row2, int col2);
1536
1537 //generate new matrix from sparse matrix, only from row1 to row2, col1 to col2 (incl.), resort with vector r
1538 void CopyFrom(const SparseMatrix& mati, int row1, int col1, int row2, int col2, const IVector& r);
1539
Init()1540 void Init()
1541 {
1542 rows=0;
1543 cols=0;
1544 mat=NULL;
1545 rowind = NULL;
1546 rowlen = NULL;
1547 colind = NULL;
1548 lalloc = 0;
1549 dnull = 0.;
1550 };
1551
1552 void GetMatrix(Matrix& m) const;
1553
GetMatrix()1554 Matrix GetMatrix() const
1555 {
1556 //very slow, only for printing!!!!!
1557 Matrix m;
1558 GetMatrix(m);
1559 return m;
1560 }
1561
1562 void SetSize(int rowsi, int colsi, int minrowsize = 1)
1563 {
1564 int resize = 0;
1565 if (rows!=rowsi || cols!=colsi) resize = 1;
1566 else
1567 {
1568 int i=0;
1569 while (!resize && i<rows)
1570 {
1571 if (RowLenMax0(i) < minrowsize) resize = 1;
1572 i++;
1573 }
1574 }
1575
1576 if (resize)
1577 {
1578 if (rowsi == 0 || colsi ==0)
1579 {
1580 rows = rowsi;
1581 cols = colsi;
1582 return;
1583 }
1584
1585 int rowsize = lalloc/rowsi; //round off!
1586 if (rowsize < minrowsize) rowsize = minrowsize;
1587
1588 int neededmem = rowsize*rowsi;
1589
1590 if (neededmem > lalloc || rowsi > rows)
1591 {
1592 Destroy();
1593
1594 lalloc = neededmem;
1595 GenMat(lalloc, rowsi);
1596 }
1597 rows = rowsi;
1598 cols = colsi;
1599
1600 int cnt = 0;
1601 for (int i = 0; i < rows; i++)
1602 {
1603 rowind[i] = cnt;
1604 cnt += rowsize;
1605 rowlen[i] = 0;
1606 }
1607 }
1608 }
1609
SetSizePerColumn(int rowsi,int colsi,const TArray<int> & size_per_column)1610 void SetSizePerColumn(int rowsi, int colsi, const TArray<int>& size_per_column)
1611 {
1612 //first it could be checked if resize is necessary!
1613 if (rowsi == 0 || colsi ==0)
1614 {
1615 rows = rowsi;
1616 cols = colsi;
1617 return;
1618 }
1619
1620 int neededmem = 0;
1621 for (int i=1; i<=size_per_column.Length(); i++)
1622 {
1623 neededmem += size_per_column(i);
1624 }
1625
1626 if (neededmem > lalloc || rowsi > rows)
1627 {
1628 Destroy();
1629
1630 lalloc = neededmem;
1631 GenMat(lalloc, rowsi);
1632 }
1633 rows = rowsi;
1634 cols = colsi;
1635
1636 int cnt = 0;
1637 for (int i = 0; i < rows; i++)
1638 {
1639 rowind[i] = cnt;
1640 cnt += size_per_column(i+1);
1641 rowlen[i] = 0;
1642 }
1643 }
1644
GetMaxlenVector()1645 Vector GetMaxlenVector() const
1646 {
1647 Vector v(rows);
1648 for (int i=0; i < rows; i++)
1649 v.Elem0(i) = RowLenMax0(i);
1650 return v;
1651 }
1652
CheckMatrix()1653 int CheckMatrix() const
1654 {
1655 int n = 0;
1656 for (int i = 0; i < rows; i++)
1657 {
1658 for (int j = 1; j < rowlen[i]; j++)
1659 {
1660 if (colind[rowind[i]+j] <= colind[rowind[i]+j-1] ||
1661 (colind[rowind[i]+j] > cols-1)) n++;
1662 }
1663 }
1664 return n;
1665 }
1666
CountZeroEntries()1667 int CountZeroEntries() const
1668 {
1669 int n = 0;
1670 for (int i = 0; i < rows; i++)
1671 {
1672 for (int j = 0; j < rowlen[i]; j++)
1673 {
1674 if (mat[rowind[i]+j] == 0) n++;
1675 }
1676 }
1677 return n;
1678 }
1679
CountEntries()1680 int CountEntries() const
1681 {
1682 int n=0;
1683 for (int i = 0; i < rows; i++)
1684 {
1685 n += rowlen[i];
1686 }
1687 return n;
1688 }
1689
1690 int EliminateZeroEntries(double mindouble=0);
1691
1692 void PrintData() const;
1693
GetMaxRowlen()1694 int GetMaxRowlen() const
1695 {
1696 int n=0;
1697 for (int i = 0; i < rows; i++)
1698 {
1699 n = Maximum(rowlen[i],n);
1700 }
1701 return n;
1702 }
1703
1704 //give each row twice space if more than 50% is used
1705 void ReAllocate();
1706
1707 //Access-functions for private data
Getrows()1708 int Getrows() const {return rows; };
Getcols()1709 int Getcols() const {return cols; };
GetLAlloc()1710 int GetLAlloc() const {return lalloc;};
1711
1712 //Assignment-operator
1713 SparseMatrix& operator= (const SparseMatrix& mati)
1714 {
1715 if (this == &mati) {return *this;}
1716
1717 if (mati.lalloc > lalloc || mati.rows > rows)
1718 {
1719 Destroy();
1720
1721 if (rows!=0 && cols!=0)
1722 {
1723 lalloc = mati.lalloc;
1724
1725 GenMat(lalloc, mati.rows);
1726 }
1727 }
1728
1729 rows=mati.rows;
1730 cols=mati.cols;
1731
1732 if (mati.mat==NULL || rows*cols==0)
1733 {
1734 mat = NULL;
1735 }
1736 else
1737 {
1738 memcpy(mat,mati.mat,rows*cols*sizeof(double));
1739 memcpy(colind,mati.colind,rows*cols*sizeof(int));
1740 memcpy(rowind,mati.rowind,rows*sizeof(int));
1741 memcpy(rowlen,mati.rowlen,rows*sizeof(int));
1742 }
1743 return *this;
1744 }
1745
1746 //fill with zeros, keep structure
FillWithZeros()1747 void FillWithZeros()
1748 {
1749 for (int i=0; i < rows; i++)
1750 {
1751 rowlen[i] = 0;
1752 }
1753 }
1754
1755 //Referencing access-operator on element using row- and column-values
operator()1756 double& operator()(int row, int col)
1757 {
1758 #ifndef __QUICKMATH
1759 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
1760 #endif
1761 //increase when necessary???
1762 int i = 0;
1763 while (i < rowlen[row-1] && colind[rowind[row-1]+i] < col-1) {i++;}
1764
1765 if (i < rowlen[row-1] && colind[rowind[row-1]+i] == col-1)
1766 {
1767 return mat[rowind[row-1]+i];
1768 }
1769 else
1770 {
1771 int ind = i;
1772
1773 if (rowlen[row-1] == RowLenMax0(row-1)) ReAllocate();
1774
1775 for (i = rowlen[row-1]-1; i >= ind; i--)
1776 {
1777 colind[rowind[row-1]+i+1] = colind[rowind[row-1]+i];
1778 mat[rowind[row-1]+i+1] = mat[rowind[row-1]+i];
1779 }
1780 rowlen[row-1]++;
1781 mat[rowind[row-1]+ind] = 0;
1782 colind[rowind[row-1]+ind] = col-1;
1783 return mat[rowind[row-1]+ind];
1784 }
1785 };
1786
1787
1788 //Referencing constant access-operator on element using row- and column-values
operator()1789 const double& operator()(int row, int col) const
1790 {
1791 #ifndef __QUICKMATH
1792 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
1793 #endif
1794 for (int i=0; i < rowlen[row-1]; i++)
1795 {
1796 if (colind[rowind[row-1]+i] == col-1)
1797 return mat[rowind[row-1]+i];
1798 }
1799 //static double dnull=0;
1800 return dnull;
1801 };
1802
1803 //Referencing constant access-operator on element using row- and column-values
Get(int row,int col)1804 const double& Get(int row, int col) const
1805 {
1806 #ifndef __QUICKMATH
1807 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
1808 #endif
1809 for (int i=0; i < rowlen[row-1]; i++)
1810 {
1811 if (colind[rowind[row-1]+i] == col-1)
1812 return mat[rowind[row-1]+i];
1813 }
1814 //static double dnull=0;
1815 return dnull;
1816 };
1817
1818 //Referencing constant access-operator on element using row- and column-values
1819 //uses last sparsecolindex of same row for faster access
Get(int row,int col,int & lastsparsecol)1820 const double& Get(int row, int col, int& lastsparsecol) const
1821 {
1822 #ifndef __QUICKMATH
1823 release_assert((row>0) && (row<=rows) && (col>0) && (col<=cols));
1824 #endif
1825 int i=lastsparsecol-1;
1826 while (i < rowlen[row-1] && colind[rowind[row-1]+i] < col-1) i++;
1827
1828 lastsparsecol = i+1; //i+1 for faster incremental access;
1829 //if (lastsparsecol >1 ) lastsparsecol--;
1830
1831 if (i < rowlen[row-1] && colind[rowind[row-1]+i] == col-1) return mat[rowind[row-1]+i];
1832
1833 //static double dnull=0;
1834 return dnull;
1835 };
1836
1837 //maximum number of elements in row-1
RowLenMax0(int row)1838 int RowLenMax0(int row) const
1839 {
1840 if (row < rows-1)
1841 {
1842 return rowind[row+1]-rowind[row];
1843 }
1844 else
1845 {
1846 return lalloc-rowind[row];
1847 }
1848 }
1849
1850 //Returns the maximum-norm (largest value in matrix)
MaxNorm()1851 double MaxNorm() const
1852 {
1853 double res=0;
1854 for(int row = 0; row < rows; row++)
1855 {
1856 for(int i = 0; i < rowlen[row]; i++)
1857 {
1858 res = Maximum(fabs(mat[rowind[row]+i]),res);
1859 }
1860 }
1861 return res;
1862 }
1863
1864 //Returns the quadratic-norm
Norm2()1865 double Norm2() const
1866 {
1867 double res=0;
1868 for(int row = 0; row < rows; row++)
1869 {
1870 for(int i = 0; i < rowlen[row]; i++)
1871 {
1872 res += Sqr(mat[rowind[row]+i]);
1873 }
1874 }
1875 return sqrt(res);
1876 }
1877
FillInFact()1878 double FillInFact() const
1879 {
1880 if (rows*cols == 0) return 0;
1881 int used = 0;
1882
1883 for (int i = 0; i < rows; i++)
1884 {
1885 used += rowlen[i];
1886 }
1887 return (double)used/(double)(rows*cols);
1888 }
1889
1890 //Returns the lower and upper bandwith of matrix
1891 void GetBandWidths(int& lb, int& ub) const;
1892
GetBandWidth()1893 int GetBandWidth() const
1894 {
1895 int a,b;
1896 GetBandWidths(a,b);
1897 return 1+Maximum(a,b);
1898 }
1899
GetBandMatrix(Matrix & bm,int lu)1900 void GetBandMatrix(Matrix& bm, int lu) const
1901 {
1902 if (0) //frage: warum funktioniert die erste version nicht???
1903 {
1904 int col_band = lu*2+1;
1905 int n = rows;
1906
1907 bm.SetSize(n,col_band);
1908
1909 int ind;
1910 int lastsparseind;
1911 for (int i = 1; i <= n; i++)
1912 {
1913 lastsparseind = 1;
1914 for (int j = 1; j <= col_band; j++)
1915 {
1916 ind = i-(lu+1)+j;
1917 if (ind >= 1 && ind <= n)
1918 {
1919 bm(i,j) = Get(i, ind);
1920 //bm(i,j) = Get(i, ind, lastsparseind);
1921 }
1922 else bm(i,j) = 0;
1923 }
1924 }
1925 }
1926 else
1927 {
1928 int col_band = lu*2+1;
1929 int n = rows;
1930
1931 bm.SetSize(n,col_band);
1932 bm.FillWithZeros();
1933
1934 int ind, k, row;
1935 for (int i = 1; i <= n; i++)
1936 {
1937 row = i-1;
1938 int rind = rowind[row];
1939 for (int j = 0; j < rowlen[row]; j++)
1940 {
1941 k = colind[rind+j]+1;
1942 ind = k-i+(lu+1);
1943 if (k >= 1 && k <= n && ind >= 1 && ind <= col_band)
1944 {
1945 bm(i,ind) = mat[rind+j];
1946 }
1947 }
1948 }
1949 }
1950 }
1951
Mult(double x)1952 void Mult(double x)
1953 {
1954 for(int row = 0; row < rows; row++)
1955 {
1956 for(int i = 0; i < rowlen[row]; i++)
1957 {
1958 mat[rowind[row]+i] *= x;
1959 }
1960 }
1961 }
1962
1963 //Add (possibly smaller) Matrix m of size (lrow x lcol) at reference position
1964 void AddMatrix(const TArray<int>& rowref, const TArray<int>& colref, int lrow, int lcol, const Matrix& m);
1965
1966 //Take Submatrix sm at (sr,sr), multiply with x and add it to local matrix at (r,c), size: (nr x nc)
1967 void AddSubmatrix(const Matrix& sm, int sr, int sc, int r, int c, int nr, int nc, double x);
1968
1969 //Take Submatrix sm at (sr,sc), multiply with x and add it to local matrix at (r,c), size: (nr x nc)
1970 void AddSubmatrix(const SparseMatrix& sm, int sr, int sc, int r, int c, int nr, int nc, double x);
1971 //set column vector, only entries of v from vr1 to vr2, insert at column col
1972 void SetColVector(const Vector& v, int vr1, int vr2, int col);
1973
1974 void SetColVector(const SparseVector& v, int col);
1975
1976 SparseMatrix& operator*= (const double& val)
1977 {
1978 //global_uo << "operator *= for sparse matrices not yet tested\n";
1979 for(int row = 0; row < rows; row++)
1980 {
1981 for(int i = 0; i < rowlen[row]; i++)
1982 {
1983 mat[rowind[row]+i] *= val;
1984 }
1985 }
1986 return *this;
1987 }
1988
1989 SparseMatrix& operator+= (const SparseMatrix& m)
1990 {
1991 //global_uo << "operator += for sparse matrices not yet tested\n";
1992 AddSubmatrix(m,1,1,1,1,rows,cols,1);
1993 return *this;
1994 }
1995
1996 SparseMatrix& operator-= (const SparseMatrix& m)
1997 {
1998 AddSubmatrix(m,1,1,1,1,rows,cols,-1);
1999 return *this;
2000 }
2001
2002 /*
2003 //Arithmetic operations with 1 parameter
2004 Matrix& operator-= (const Matrix& m);
2005
2006 friend void Mult(const SparseMatrix& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2
2007 friend void Mult(const Matrix3D& m1, const SparseMatrix& m2, Matrix& res); //computes res=m1*m2
2008 */
2009
2010 friend Matrix operator* (const SparseMatrix& m1, const Matrix& m2);
2011 friend void Mult(const SparseMatrix& m, const Vector& v, Vector& res); //computes res=m*v
2012 friend void Mult(const SparseMatrix& m1, const Matrix& m2, Matrix& res); //computes res=m1*m2
2013
2014 friend void Mult(const SparseMatrix& m1, const SparseMatrix& m2, SparseMatrix& res); //computes res=m1*m2
2015 virtual void ComputeTranspose(SparseMatrix& res) const; //computes res=m1^T
2016 friend double MultRow(const SparseMatrix& m1, const SparseMatrix& m2, int row1, int row2); //mult row1 of m1 times row2 of m2
2017
2018
2019 void PrintToMatlabSparse(ostream& os);
2020
2021 //writes out matrix with constandt width and factor normalized
2022 friend ostream& operator<<(ostream& os, const SparseMatrix& m);
2023
GetMatPtrs(int * & colindI,int * & rowptrI,double * & matI)2024 void GetMatPtrs(int*& colindI, int*& rowptrI, double*& matI) const
2025 {
2026 colindI = colind;
2027 rowptrI = rowind;
2028 matI = mat;
2029 }
GetRowLenArray()2030 int* GetRowLenArray() const {return rowlen;}
2031
RowLen(int i)2032 virtual int RowLen(int i) const {return rowlen[i-1];}
GetRowEntry(int row,int i)2033 virtual double GetRowEntry(int row, int i) const {return mat[rowind[row-1]+i-1];}
GetRowColind(int row,int i)2034 virtual int GetRowColind(int row, int i) const {return colind[rowind[row-1]+i-1];}
2035
2036 private:
2037
2038 int rows, cols;
2039 int lalloc;
2040
2041 double* mat;
2042 int* colind;
2043 int* rowind;
2044 int* rowlen;
2045
2046 double dnull;
2047 };
2048
2049 //$ PG 2013-10-2: deleted class HarwellBoeingMatrix
HarwellBoeingMatrix()2050 class HarwellBoeingMatrix { public: HarwellBoeingMatrix() {assert(0 && "ERROR: HarwellBoeingMatrix removed");} };
2051
2052 //struct SuperMatrix;
2053
2054
2055 // AP
2056 class SuperLUInverse
2057 {
2058 public:
SuperLUInverse()2059 SuperLUInverse() : sparseA(0), rows(0), cols(0), perm_r(0), perm_c(0), isFirstStep(1), isFactorized(0),
2060 SuperA(0), SuperL(0), SuperU(0)
2061 {
2062 ;
2063 }
2064
SuperLUInverse(const SparseMatrix & matA)2065 SuperLUInverse:: SuperLUInverse(const SparseMatrix & matA) : sparseA(&matA), isFactorized(0), isFirstStep(1),
2066 SuperA(0), SuperL(0), SuperU(0)
2067 {
2068 sparseA =&matA;
2069 rows = matA.Getrows();
2070 cols = matA.Getcols();
2071
2072 isFactorized = 0;
2073 isFirstStep = 1;
2074 perm_r = new int(rows);
2075 perm_c = new int(cols);
2076 }
2077
2078
2079 ~SuperLUInverse();
2080
SetSparseMatrix(const SparseMatrix & matA)2081 void SetSparseMatrix(const SparseMatrix& matA)
2082 {
2083 sparseA =&matA;
2084 rows = matA.Getrows();
2085 cols = matA.Getcols();
2086
2087 isFactorized = 0;
2088 isFirstStep = 1;
2089 delete[] perm_r;
2090 delete[] perm_c;
2091 perm_r = new int(rows);
2092 perm_c = new int(cols);
2093 }
2094
2095 int Solve(Vector& q);
2096 int Factorize();
2097 int SuperLUApply(Vector& q);
Apply(Vector & q)2098 int Apply(Vector& q) { return SuperLUApply(q);}
2099
IsFirstStep()2100 int& IsFirstStep() {return isFirstStep; }
IsFirstStep()2101 const int& IsFirstStep() const {return isFirstStep; }
2102
2103
2104 private:
2105 // pointer to original sparsematrix
2106 const SparseMatrix* sparseA;
2107 // size, has to be quadratic
2108 int rows, cols;
2109 // SuperLU-Data
2110 SuperMatrix* SuperA;
2111 SuperMatrix* SuperU;
2112 SuperMatrix* SuperL;
2113 int* perm_r;
2114 int* perm_c;
2115 int isFirstStep;
2116 // is factorization done?
2117 int isFactorized;
2118
2119 };
2120 // END AP
2121
2122
2123 //solve A*q = rhs; q contains right-hand-side on input and the solution on output
2124 int SuperLUSolve(SparseMatrix& A, Vector& q);
2125 int SuperLUFactorize(SparseMatrix& A, SuperMatrix*& SuperA, SuperMatrix*& L, SuperMatrix*& U, int*& perm_r, int*& perm_c, int isFirstStep = 1);
2126 int SuperLUApply(SparseMatrix& A, SuperMatrix*& SuperA, SuperMatrix*& L, SuperMatrix*& U, int*& perm_r, int*& perm_c, Vector& q);
2127
2128
2129 //$ DR 2013-10-04:[ the following functions are not available anymore
2130 //int banddec(int n, int ld, int ud, Matrix& pmat, IVector& perm, int& signd);
2131 //int bandsol(int n, int ld, int ud, Matrix& pmat, Vector& b, const IVector& perm);
2132 //int Band_LU_Dec(Matrix& a, int n, int m1, int m2, Matrix& al, IVector& indx, double& d);
2133 //void Band_LU_Bks(const Matrix& a, int n, int m1, int m2, const Matrix& al, const IVector& indx, Vector& b);
2134 //void Band_LU_Bks0(const Matrix& a, int n, int m1, int m2, const Matrix& al, const IVector& indx, Vector& b);
2135 //$ DR 2013-10-04:]
2136
2137 class PardisoInverse;
2138
2139 // The SparseInverse uses Pardiso or SuperLU internally
2140 class SparseInverse
2141 {
2142 public:
SparseInverse()2143 SparseInverse() : pardisoinv(NULL), superluinv(NULL)
2144 {
2145 }
2146
2147 SparseInverse(const SparseMatrix & matA) ;
2148
2149 ~SparseInverse();
2150 void SetSparseMatrix(const SparseMatrix& matA);
2151 int Solve(Vector& q);
2152 int Factorize();
2153 int Apply(Vector& q);
2154 int& IsFirstStep() ;
2155 const int& IsFirstStep() const ;
2156
2157 private:
2158 SuperLUInverse* superluinv;
2159 PardisoInverse* pardisoinv;
2160 };
2161
2162
2163 //symmetrical smoothing
2164 void SmoothenDoubleArray(Vector& values, Vector& result, int smoothing_length);
2165
2166
2167 // this function computes the local maxima and minima values if a vector 'data' from left to right (in case of similar points, the first point from left to right is used - see '�' below)
2168 // if only max-values are computed, set flag compute_minima to zero and use dummys for TArrays 'min_indices', and 'min_values'
2169 //
2170 //max1
2171 //o max2
2172 // \ �-- maxN
2173 // \ / \ max3 �
2174 // � � � .../
2175 // \ / \ /
2176 // � �
2177 // min1 min2
2178 void GetExtremeValues(const Vector& data, TArray<int>& max_indices,TArray<double>& max_values, TArray<int>& min_indices,TArray<double>& min_values, int compute_minima=1);