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