Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
UblasCustomFunctions.hpp
Go to the documentation of this file.
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifndef UBLASCUSTOMFUNCTIONS_HPP_
37#define UBLASCUSTOMFUNCTIONS_HPP_
38
45#include "UblasIncludes.hpp"
46
47#include "Exception.hpp"
49
50// COMMON DETERMINANTS - SQUARE
51
58template<class T>
59inline T Determinant(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM)
60{
61 using namespace boost::numeric::ublas;
62
63 return rM(0,0);
64}
65
72template<class T>
73T Determinant(const boost::numeric::ublas::c_matrix<T,2,2>& rM)
74{
75 using namespace boost::numeric::ublas;
76
77 return rM(0,0)*rM(1,1) - rM(1,0)*rM(0,1);
78}
79
86template<class T>
87T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM)
88{
89 using namespace boost::numeric::ublas;
90
91 return rM(0,0) * (rM(1,1)*rM(2,2) - rM(1,2)*rM(2,1))
92 - rM(0,1) * (rM(1,0)*rM(2,2) - rM(1,2)*rM(2,0))
93 + rM(0,2) * (rM(1,0)*rM(2,1) - rM(1,1)*rM(2,0));
94}
95
96// COMMON GENERALIZED DETERMINANTS - NOT SQUARE
97
105template<class T>
106T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM)
107{
108 using namespace boost::numeric::ublas;
109 c_matrix<T,2,2> product = prod(trans(rM), rM);
110 return std::sqrt(Determinant(product));
111}
112
120template<class T>
121T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM)
122{
123 using namespace boost::numeric::ublas;
124 return std::sqrt(rM(0,0)*rM(0,0) + rM(1,0)*rM(1,0) + rM(2,0)*rM(2,0));
125}
126
134template<class T>
135T Determinant(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM)
136{
137 using namespace boost::numeric::ublas;
138 return std::sqrt(rM(0,0) * rM(0,0) + rM(1,0) * rM(1,0));
139}
140
147template<class T>
148T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 0>& rM)
149{
151}
152
159template<class T>
160T Determinant(const boost::numeric::ublas::c_matrix<T, 2, 0>& rM)
161{
163}
164
171template<class T>
172T Determinant(const boost::numeric::ublas::c_matrix<T, 1, 0>& rM)
173{
175}
176
177// COMMON SUBDETERMINANTS - SQUARE
178
188template<class T>
189T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM, const unsigned missrow, const unsigned misscol)
190{
191 using namespace boost::numeric::ublas;
192
193 assert(missrow == 0);
194 assert(misscol == 0);
195 return 1.0;
196}
197
206template<class T>
207T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 2>& rM, const unsigned missrow, const unsigned misscol)
208{
209 using namespace boost::numeric::ublas;
210
211 assert(missrow < 2);
212 assert(misscol < 2);
213
214 unsigned row = (missrow==1) ? 0 : 1;
215 unsigned col = (misscol==1) ? 0 : 1;
216 return rM(row,col);
217}
218
227template<class T>
228T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM, const unsigned missrow, const unsigned misscol)
229{
230 using namespace boost::numeric::ublas;
231
232 assert(missrow < 3);
233 assert(misscol < 3);
234
235 unsigned lorow = (missrow==0) ? 1 : 0;
236 unsigned hirow = (missrow==2) ? 1 : 2;
237 unsigned locol = (misscol==0) ? 1 : 0;
238 unsigned hicol = (misscol==2) ? 1 : 2;
239 return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
240}
241
242// COMMON SUBDETERMINANTS - NOT SQUARE
243
252template<class T>
253T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM, const unsigned missrow, const unsigned misscol)
254{
255 using namespace boost::numeric::ublas;
256
257 assert(missrow < 3);
258 //assert(misscol < 2); //Don't assert this since it is used
259
260 unsigned lorow = (missrow==0) ? 1 : 0;
261 unsigned hirow = (missrow==2) ? 1 : 2;
262 unsigned locol = (misscol==0) ? 1 : 0;
263 unsigned hicol = (misscol==2) ? 1 : 2;
264 return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
265}
266
275template<class T>
276T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM, const unsigned missrow, const unsigned misscol)
277{
278 using namespace boost::numeric::ublas;
279
280 assert(missrow < 3);
281 assert(misscol < 1);
282
283 unsigned lorow = (missrow==0) ? 1 : 0;
284 unsigned hirow = (missrow==2) ? 1 : 2;
285 unsigned locol = (misscol==0) ? 1 : 0;
286 unsigned hicol = (misscol==2) ? 1 : 2;
287 return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
288}
289
298template<class T>
299T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM, const unsigned missrow, const unsigned misscol)
300{
301 using namespace boost::numeric::ublas;
302
303 assert(missrow < 2);
304 assert(misscol < 1);
305
306 unsigned row = (missrow==1) ? 0 : 1;
307 unsigned col = (misscol==1) ? 0 : 1;
308 return rM(row,col);
309}
310
311#if defined(__xlC__)
312/* IBM compiler doesn't support zero-sized arrays*/
313#else //#if defined(__xlC__)
322template<class T>
323T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 0>& rM, const unsigned missrow, const unsigned misscol)
324{
326}
327
336template<class T>
337T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 0>& rM, const unsigned missrow, const unsigned misscol)
338{
340}
341
350template<class T>
351T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 1, 0>& rM, const unsigned missrow, const unsigned misscol)
352{
354}
355#endif //#if defined(__xlC__)
356
357// COMMON INVERSES - SQUARE
358
366template<class T>
367boost::numeric::ublas::c_matrix<T, 1, 1> Inverse(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM)
368{
369 using namespace boost::numeric::ublas;
370
371 c_matrix<T,1,1> inverse;
372 T det = Determinant(rM);
373 assert(fabs(det) > DBL_EPSILON); // else it is a singular matrix
374 inverse(0,0) = 1.0/det;
375 return inverse;
376}
377
385template<class T>
386boost::numeric::ublas::c_matrix<T, 2, 2> Inverse(const boost::numeric::ublas::c_matrix<T, 2, 2>& rM)
387{
388 using namespace boost::numeric::ublas;
389
390 c_matrix<T, 2, 2> inverse;
391 T det = Determinant(rM);
392
393 assert( fabs(det) > DBL_EPSILON ); // else it is a singular matrix
394 inverse(0,0) = rM(1,1)/det;
395 inverse(0,1) = -rM(0,1)/det;
396 inverse(1,0) = -rM(1,0)/det;
397 inverse(1,1) = rM(0,0)/det;
398 return inverse;
399}
400
408template<class T>
409boost::numeric::ublas::c_matrix<T, 3, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM)
410{
411 using namespace boost::numeric::ublas;
412
413 c_matrix<T, 3, 3> inverse;
414 T det = Determinant(rM);
415 assert(fabs(det) > DBL_EPSILON); // else it is a singular matrix
416
417 inverse(0,0) = (rM(1,1)*rM(2,2) - rM(1,2)*rM(2,1))/det;
418 inverse(1,0) = -(rM(1,0)*rM(2,2) - rM(1,2)*rM(2,0))/det;
419 inverse(2,0) = (rM(1,0)*rM(2,1) - rM(1,1)*rM(2,0))/det;
420 inverse(0,1) = -(rM(0,1)*rM(2,2) - rM(0,2)*rM(2,1))/det;
421 inverse(1,1) = (rM(0,0)*rM(2,2) - rM(0,2)*rM(2,0))/det;
422 inverse(2,1) = -(rM(0,0)*rM(2,1) - rM(0,1)*rM(2,0))/det;
423 inverse(0,2) = (rM(0,1)*rM(1,2) - rM(0,2)*rM(1,1))/det;
424 inverse(1,2) = -(rM(0,0)*rM(1,2) - rM(0,2)*rM(1,0))/det;
425 inverse(2,2) = (rM(0,0)*rM(1,1) - rM(0,1)*rM(1,0))/det;
426
427 return inverse;
428}
429
430// COMMON PSEUDO-INVERSES - NOT SQUARE
431
439template<class T>
440boost::numeric::ublas::c_matrix<T, 2, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM)
441{
442 using namespace boost::numeric::ublas;
443
444 c_matrix<T, 2, 3> inverse;
445
446 //
447 // calculate (T'T)^-1, where T'T = (a b)
448 // (c d)
449
450 T a = rM(0,0)*rM(0,0) + rM(1,0)*rM(1,0) + rM(2,0)*rM(2,0);
451 T b = rM(0,0)*rM(0,1) + rM(1,0)*rM(1,1) + rM(2,0)*rM(2,1);
452 T c = b;
453 T d = rM(0,1)*rM(0,1) + rM(1,1)*rM(1,1) + rM(2,1)*rM(2,1);
454
455 T det = a*d - b*c;
456
457 T a_inv = d/det;
458 T b_inv = -b/det;
459 T c_inv = -c/det;
460 T d_inv = a/det;
461
462 inverse(0,0) = a_inv*rM(0,0) + b_inv*rM(0,1);
463 inverse(1,0) = c_inv*rM(0,0) + d_inv*rM(0,1);
464 inverse(0,1) = a_inv*rM(1,0) + b_inv*rM(1,1);
465 inverse(1,1) = c_inv*rM(1,0) + d_inv*rM(1,1);
466 inverse(0,2) = a_inv*rM(2,0) + b_inv*rM(2,1);
467 inverse(1,2) = c_inv*rM(2,0) + d_inv*rM(2,1);
468
469 return inverse;
470}
471
479template<class T>
480boost::numeric::ublas::c_matrix<T, 1, 2> Inverse(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM)
481{
482 using namespace boost::numeric::ublas;
483
484 c_matrix<T, 1, 2> inverse;
485 T det = Determinant(rM);
486
487 inverse(0,0) = rM(0,0)/det/det;
488 inverse(0,1) = rM(1,0)/det/det;
489
490 return inverse;
491}
492
500template<class T>
501boost::numeric::ublas::c_matrix<T, 1, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM)
502{
503 using namespace boost::numeric::ublas;
504
505 c_matrix<T, 1, 3> inverse;
506 T det = Determinant(rM);
507
508 inverse(0,0) = rM(0,0)/det/det;
509 inverse(0,1) = rM(1,0)/det/det;
510 inverse(0,2) = rM(2,0)/det/det;
511
512 return inverse;
513}
514
515// COMMON MATRIX TRACES
516
523template<class T>
524inline T Trace(const c_matrix<T, 1, 1>& rM)
525{
526 return rM(0,0);
527}
528
535template<class T>
536inline T Trace(const c_matrix<T, 2, 2>& rM)
537{
538 return rM(0,0) + rM(1,1);
539}
540
547template<class T>
548inline T Trace(const c_matrix<T, 3, 3>& rM)
549{
550 return rM(0,0) + rM(1,1) + rM(2,2);
551}
552
559template<class T>
560inline T Trace(const c_matrix<T, 4, 4>& rM)
561{
562 return rM(0,0) + rM(1,1) + rM(2,2) + rM(3,3);
563}
564
565// OTHER MATRIX FUNCTIONS (INVARIANTS, EIGENVECTORS)
566
575template<class T>
576T SecondInvariant(const c_matrix<T, 3, 3>& rM)
577{
578 return rM(0,0)*rM(1,1) + rM(1,1)*rM(2,2) + rM(2,2)*rM(0,0)
579 - rM(1,0)*rM(1,0) - rM(2,1)*rM(2,1) - rM(2,0)*rM(2,0);
580}
581
589template<class T>
590T SecondInvariant(const c_matrix<T, 2, 2>& rM)
591{
592 return Determinant(rM);
593}
594
604c_vector<double,3> CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix<double, 3, 3>& rA);
605
612double CalculateMaxEigenpair(c_matrix<double, 3, 3>& rA, c_vector<double, 3>& rEigenvector);
613
614 //COMMON VECTOR FUNCTIONS
615
616
625template<class T>
626c_vector<T, 1> VectorProduct(const c_vector<T, 1>& rA, const c_vector<T, 1>& rB)
627{
629}
638template<class T>
639c_vector<T, 2> VectorProduct(const c_vector<T, 2>& rA, const c_vector<T, 2>& rB)
640{
642}
643
651template<class T>
652c_vector<T, 3> VectorProduct(const c_vector<T, 3>& rA, const c_vector<T, 3>& rB)
653{
654
655 c_vector<T, 3> result;
656
657 result(0) = rA(1)*rB(2) - rA(2)*rB(1);
658 result(1) = rA(2)*rB(0) - rA(0)*rB(2);
659 result(2) = rA(0)*rB(1) - rA(1)*rB(0);
660
661 return result;
662}
663
670c_vector<double, 1> Create_c_vector(double x);
671
679c_vector<double, 2> Create_c_vector(double x, double y);
680
689c_vector<double, 3> Create_c_vector(double x, double y, double z);
690
691#endif /*UBLASCUSTOMFUNCTIONS_HPP_*/
#define NEVER_REACHED
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
T SecondInvariant(const c_matrix< T, 3, 3 > &rM)
T SubDeterminant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM, const unsigned missrow, const unsigned misscol)
T Trace(const c_matrix< T, 1, 1 > &rM)
c_vector< double, 3 > CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix< double, 3, 3 > &rA)
c_vector< double, 1 > Create_c_vector(double x)
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
double CalculateMaxEigenpair(c_matrix< double, 3, 3 > &rA, c_vector< double, 3 > &rEigenvector)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)