2.5: Remove OOPS code from the outliner space, as discussed
[blender-staging.git] / extern / ode / dist / ode / fbuild / test_ldlt.cpp
1 /*************************************************************************
2  *                                                                       *
3  * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith.       *
4  * All rights reserved.  Email: russ@q12.org   Web: www.q12.org          *
5  *                                                                       *
6  * This library is free software; you can redistribute it and/or         *
7  * modify it under the terms of EITHER:                                  *
8  *   (1) The GNU Lesser General Public License as published by the Free  *
9  *       Software Foundation; either version 2.1 of the License, or (at  *
10  *       your option) any later version. The text of the GNU Lesser      *
11  *       General Public License is included with this library in the     *
12  *       file LICENSE.TXT.                                               *
13  *   (2) The BSD-style license that is included with this library in     *
14  *       the file LICENSE-BSD.TXT.                                       *
15  *                                                                       *
16  * This library is distributed in the hope that it will be useful,       *
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files    *
19  * LICENSE.TXT and LICENSE-BSD.TXT for more details.                     *
20  *                                                                       *
21  *************************************************************************/
22
23 #include <stdio.h>
24 #include <malloc.h>
25 #include "ode/ode.h"
26
27 #define ALLOCA dALLOCA16
28
29 //****************************************************************************
30 // constants
31
32 #ifdef dSINGLE
33 #define TOL (1e-4)
34 #else
35 #define TOL (1e-10)
36 #endif
37
38 //****************************************************************************
39 // test L*X=B solver accuracy.
40
41 void testSolverAccuracy (int n)
42 {
43   int i;
44   int npad = dPAD(n);
45   dReal *L = (dReal*) ALLOCA (n*npad*sizeof(dReal));
46   dReal *B = (dReal*) ALLOCA (n*sizeof(dReal));
47   dReal *B2 = (dReal*) ALLOCA (n*sizeof(dReal));
48   dReal *X = (dReal*) ALLOCA (n*sizeof(dReal));
49
50   // L is a random lower triangular matrix with 1's on the diagonal
51   dMakeRandomMatrix (L,n,n,1.0);
52   dClearUpperTriangle (L,n);
53   for (i=0; i<n; i++) L[i*npad+i] = 1;
54
55   // B is the right hand side
56   dMakeRandomMatrix (B,n,1,1.0);
57   memcpy (X,B,n*sizeof(dReal)); // copy B to X
58
59   dSolveL1 (L,X,n,npad);
60
61   /*
62   dPrintMatrix (L,n,n);
63   printf ("\n");
64   dPrintMatrix (B,n,1);
65   printf ("\n");
66   dPrintMatrix (X,n,1);
67   printf ("\n");
68   */
69
70   dSetZero (B2,n);
71   dMultiply0 (B2,L,X,n,n,1);
72   dReal error = dMaxDifference (B,B2,1,n);
73   if (error > TOL) {
74     printf ("error = %e, size = %d\n",error,n);
75   }
76 }
77
78 //****************************************************************************
79 // test L^T*X=B solver accuracy.
80
81 void testTransposeSolverAccuracy (int n)
82 {
83   int i;
84   int npad = dPAD(n);
85   dReal *L = (dReal*) ALLOCA (n*npad*sizeof(dReal));
86   dReal *B = (dReal*) ALLOCA (n*sizeof(dReal));
87   dReal *B2 = (dReal*) ALLOCA (n*sizeof(dReal));
88   dReal *X = (dReal*) ALLOCA (n*sizeof(dReal));
89
90   // L is a random lower triangular matrix with 1's on the diagonal
91   dMakeRandomMatrix (L,n,n,1.0);
92   dClearUpperTriangle (L,n);
93   for (i=0; i<n; i++) L[i*npad+i] = 1;
94
95   // B is the right hand side
96   dMakeRandomMatrix (B,n,1,1.0);
97   memcpy (X,B,n*sizeof(dReal)); // copy B to X
98
99   dSolveL1T (L,X,n,npad);
100
101   dSetZero (B2,n);
102   dMultiply1 (B2,L,X,n,n,1);
103   dReal error = dMaxDifference (B,B2,1,n);
104   if (error > TOL) {
105     printf ("error = %e, size = %d\n",error,n);
106   }
107 }
108
109 //****************************************************************************
110 // test L*D*L' factorizer accuracy.
111
112 void testLDLTAccuracy (int n)
113 {
114   int i,j;
115   int npad = dPAD(n);
116   dReal *A = (dReal*) ALLOCA (n*npad*sizeof(dReal));
117   dReal *L = (dReal*) ALLOCA (n*npad*sizeof(dReal));
118   dReal *d = (dReal*) ALLOCA (n*sizeof(dReal));
119   dReal *Atest = (dReal*) ALLOCA (n*npad*sizeof(dReal));
120   dReal *DL = (dReal*) ALLOCA (n*npad*sizeof(dReal));
121
122   dMakeRandomMatrix (A,n,n,1.0);
123   dMultiply2 (L,A,A,n,n,n);
124   memcpy (A,L,n*npad*sizeof(dReal));
125   dSetZero (d,n);
126
127   dFactorLDLT (L,d,n,npad);
128
129   // make L lower triangular, and convert d into diagonal of D
130   dClearUpperTriangle (L,n);
131   for (i=0; i<n; i++) L[i*npad+i] = 1;
132   for (i=0; i<n; i++) d[i] = 1.0/d[i];
133
134   // form Atest = L*D*L'
135   dSetZero (Atest,n*npad);
136   dSetZero (DL,n*npad);
137   for (i=0; i<n; i++) {
138     for (j=0; j<n; j++) DL[i*npad+j] = L[i*npad+j] * d[j];
139   }
140   dMultiply2 (Atest,L,DL,n,n,n);
141   dReal error = dMaxDifference (A,Atest,n,n);
142   if (error > TOL) {
143     printf ("error = %e, size = %d\n",error,n);
144   }
145
146   /*
147   printf ("\n");
148   dPrintMatrix (A,n,n);
149   printf ("\n");
150   dPrintMatrix (L,n,n);
151   printf ("\n");
152   dPrintMatrix (d,1,n);
153   */
154 }
155
156 //****************************************************************************
157 // test L*D*L' factorizer speed.
158
159 void testLDLTSpeed (int n)
160 {
161   int npad = dPAD(n);
162
163   // allocate A
164   dReal *A = (dReal*) ALLOCA (n*npad*sizeof(dReal));
165
166   // make B a symmetric positive definite matrix
167   dMakeRandomMatrix (A,n,n,1.0);
168   dReal *B = (dReal*) ALLOCA (n*npad*sizeof(dReal));
169   dSetZero (B,n*npad);
170   dMultiply2 (B,A,A,n,n,n);
171
172   // make d
173   dReal *d = (dReal*) ALLOCA (n*sizeof(dReal));
174   dSetZero (d,n);
175
176   // time several factorizations, return the minimum timing
177   double mintime = 1e100;
178   dStopwatch sw;
179   for (int i=0; i<100; i++) {
180     memcpy (A,B,n*npad*sizeof(dReal));
181     dStopwatchReset (&sw);
182     dStopwatchStart (&sw);
183
184     dFactorLDLT (A,d,n,npad);
185
186     dStopwatchStop (&sw);
187     double time = dStopwatchTime (&sw);
188     if (time < mintime) mintime = time;
189   }
190
191   printf ("%.0f",mintime * dTimerTicksPerSecond());
192 }
193
194 //****************************************************************************
195 // test solver speed.
196
197 void testSolverSpeed (int n, int transpose)
198 {
199   int i;
200   int npad = dPAD(n);
201
202   // allocate L,B,X
203   dReal *L = (dReal*) ALLOCA (n*npad*sizeof(dReal));
204   dReal *B = (dReal*) ALLOCA (n*sizeof(dReal));
205   dReal *X = (dReal*) ALLOCA (n*sizeof(dReal));
206
207   // L is a random lower triangular matrix with 1's on the diagonal
208   dMakeRandomMatrix (L,n,n,1.0);
209   dClearUpperTriangle (L,n);
210   for (i=0; i<n; i++) L[i*npad+i] = 1;
211
212   // B is the right hand side
213   dMakeRandomMatrix (B,n,1,1.0);
214
215   // time several factorizations, return the minimum timing
216   double mintime = 1e100;
217   dStopwatch sw;
218   for (int i=0; i<100; i++) {
219     memcpy (X,B,n*sizeof(dReal));       // copy B to X
220
221     if (transpose) {
222       dStopwatchReset (&sw);
223       dStopwatchStart (&sw);
224       dSolveL1T (L,X,n,npad);
225       dStopwatchStop (&sw);
226     }
227     else {
228       dStopwatchReset (&sw);
229       dStopwatchStart (&sw);
230       dSolveL1 (L,X,n,npad);
231       dStopwatchStop (&sw);
232     }
233
234     double time = dStopwatchTime (&sw);
235     if (time < mintime) mintime = time;
236   }
237
238   printf ("%.0f",mintime * dTimerTicksPerSecond());
239 }
240
241 //****************************************************************************
242 // the single command line argument is 'f' to test and time the factorizer,
243 // or 's' to test and time the solver.
244
245
246 void testAccuracy (int n, char type)
247 {
248   if (type == 'f') testLDLTAccuracy (n);
249   if (type == 's') testSolverAccuracy (n);
250   if (type == 't') testTransposeSolverAccuracy (n);
251 }
252
253
254 void testSpeed (int n, char type)
255 {
256   if (type == 'f') testLDLTSpeed (n);
257   if (type == 's') testSolverSpeed (n,0);
258   if (type == 't') testSolverSpeed (n,1);
259 }
260
261
262 int main (int argc, char **argv)
263 {
264   if (argc != 2 || argv[1][0] == 0 || argv[1][1] != 0 ||
265       (argv[1][0] != 'f' && argv[1][0] != 's' && argv[1][0] != 't')) {
266     fprintf (stderr,"Usage: test_ldlt [f|s|t]\n");
267     exit (1);
268   }
269   char type = argv[1][0];
270
271   // accuracy test: test all sizes up to 20 then all prime sizes up to 101
272   int i;
273   for (i=1; i<20; i++) {
274     testAccuracy (i,type);
275   }
276   testAccuracy (23,type);
277   testAccuracy (29,type);
278   testAccuracy (31,type);
279   testAccuracy (37,type);
280   testAccuracy (41,type);
281   testAccuracy (43,type);
282   testAccuracy (47,type);
283   testAccuracy (53,type);
284   testAccuracy (59,type);
285   testAccuracy (61,type);
286   testAccuracy (67,type);
287   testAccuracy (71,type);
288   testAccuracy (73,type);
289   testAccuracy (79,type);
290   testAccuracy (83,type);
291   testAccuracy (89,type);
292   testAccuracy (97,type);
293   testAccuracy (101,type);
294
295   // test speed on a 127x127 matrix
296   testSpeed (127,type);
297
298   return 0;
299 }