doxygen: intern/smoke tagged.
[blender.git] / intern / smoke / intern / LU_HELPER.cpp
1 /** \file smoke/intern/LU_HELPER.cpp
2  *  \ingroup smoke
3  */
4
5 #include "LU_HELPER.h"
6
7 int isNonsingular (sLU LU_) {
8       for (int j = 0; j < 3; j++) {
9          if (LU_.values[j][j] == 0)
10             return 0;
11       }
12       return 1;
13 }
14
15 sLU computeLU( float a[3][3])
16 {
17         sLU result;
18         int m=3;
19         int n=3;
20
21         //float LU_[3][3];
22         for (int i = 0; i < m; i++) {
23                         result.piv[i] = i;
24                         for (int j = 0; j < n; j++) result.values[i][j]=a[i][j];
25       }
26
27       result.pivsign = 1;
28       //Real *LUrowi = 0;;
29       //Array1D<Real> LUcolj(m);
30           //float *LUrowi = 0;
31           float LUcolj[3];
32
33       // Outer loop.
34
35       for (int j = 0; j < n; j++) {
36
37          // Make a copy of the j-th column to localize references.
38
39          for (int i = 0; i < m; i++) {
40             LUcolj[i] = result.values[i][j];
41          }
42
43          // Apply previous transformations.
44
45          for (int i = 0; i < m; i++) {
46                          //float LUrowi[3];
47                          //LUrowi = result.values[i];
48
49             // Most of the time is spent in the following dot product.
50
51             int kmax = min(i,j);
52             double s = 0.0;
53             for (int k = 0; k < kmax; k++) {
54                s += result.values[i][k]*LUcolj[k];
55             }
56
57             result.values[i][j] = LUcolj[i] -= s;
58          }
59    
60          // Find pivot and exchange if necessary.
61
62          int p = j;
63          for (int i = j+1; i < m; i++) {
64             if (abs(LUcolj[i]) > abs(LUcolj[p])) {
65                p = i;
66             }
67          }
68          if (p != j) {
69                     int k=0;
70             for (k = 0; k < n; k++) {
71                double t = result.values[p][k]; 
72                            result.values[p][k] = result.values[j][k]; 
73                            result.values[j][k] = t;
74             }
75             k = result.piv[p]; 
76                         result.piv[p] = result.piv[j]; 
77                         result.piv[j] = k;
78             result.pivsign = -result.pivsign;
79          }
80
81          // Compute multipliers.
82          
83          if ((j < m) && (result.values[j][j] != 0.0)) {
84             for (int i = j+1; i < m; i++) {
85                result.values[i][j] /= result.values[j][j];
86             }
87          }
88       }
89
90           return result;
91 }
92
93 void solveLU3x3(sLU& A, float x[3], float b[3])
94 {
95   //TNT::Array1D<float> jamaB = TNT::Array1D<float>(3, &b[0]);
96   //TNT::Array1D<float> jamaX = A.solve(jamaB);
97
98         
99   // Solve A, B
100
101         {
102       if (!isNonsingular(A)) {
103         x[0]=0.0f;
104                 x[1]=0.0f;
105                 x[2]=0.0f;
106                 return;
107       }
108
109
110           //Array1D<Real> Ax = permute_copy(b, piv);
111           float Ax[3];
112
113     // permute copy: b , A.piv
114         {
115          for (int i = 0; i < 3; i++) 
116                Ax[i] = b[A.piv[i]];
117         }
118
119       // Solve L*Y = B(piv)
120       for (int k = 0; k < 3; k++) {
121          for (int i = k+1; i < 3; i++) {
122                Ax[i] -= Ax[k]*A.values[i][k];
123             }
124          }
125       
126           // Solve U*X = Y;
127       for (int k = 2; k >= 0; k--) {
128             Ax[k] /= A.values[k][k];
129                 for (int i = 0; i < k; i++) 
130                 Ax[i] -= Ax[k]*A.values[i][k];
131       }
132      
133
134                 x[0] = Ax[0];
135                 x[1] = Ax[1];
136                 x[2] = Ax[2];
137       return;
138         }
139 }