fd74faf18df8b9e7ffc0f01a7ebf46268b6dd3cc
[blender-staging.git] / extern / ode / dist / ode / fbuild / BuildLDLT
1 #!/usr/bin/perl
2 #
3 #########################################################################
4 #                                                                       #
5 # Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith.       #
6 # All rights reserved.  Email: russ@q12.org   Web: www.q12.org          #
7 #                                                                       #
8 # This library is free software; you can redistribute it and/or         #
9 # modify it under the terms of EITHER:                                  #
10 #   (1) The GNU Lesser General Public License as published by the Free  #
11 #       Software Foundation; either version 2.1 of the License, or (at  #
12 #       your option) any later version. The text of the GNU Lesser      #
13 #       General Public License is included with this library in the     #
14 #       file LICENSE.TXT.                                               #
15 #   (2) The BSD-style license that is included with this library in     #
16 #       the file LICENSE-BSD.TXT.                                       #
17 #                                                                       #
18 # This library is distributed in the hope that it will be useful,       #
19 # but WITHOUT ANY WARRANTY; without even the implied warranty of        #
20 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files    #
21 # LICENSE.TXT and LICENSE-BSD.TXT for more details.                     #
22 #                                                                       #
23 #########################################################################
24
25 #
26 # triangular matrix solver and factorizer code generator.
27 #
28 # SOLVER
29 # ------
30 #
31 # if L is an n x n lower triangular matrix (with ones on the diagonal), the
32 # solver solves L*X=B where X and B are n x m matrices. this is the core
33 # step in L*D*L' factorization. the algorithm is (in matlab):
34 #
35 #     for i=1:n
36 #       for j=1:m
37 #         X(i,j) = B(i,j) - L(i,1:i-1)*X(1:i-1,j);
38 #       end
39 #     end
40 #
41 # note that the ordering of the (i,j) loop is somewhat arbitrary. the only
42 # prerequisite to calculating element (i,j) of X is that all X(1:i-1,j) have
43 # have already been calcuated. this gives us some flexibility.
44 #
45 # the code generated below calculates X in N1 x N1 blocks. to speed up the
46 # innermost dot product loop, the outer product trick is used. for instance,
47 # to calculate the value of the 2x2 matrix ABCD below we first iterate over
48 # the vectors (a,b,c,d) and (e,f,g,h), computing ABCD = a*e+b*f+c*g+d*h.
49 # then A and B contain the dot product values needed in the algorithm, and
50 # C and D have almost all of it. the outer product trick reduces the number
51 # of memory loads required. in this example 16 loads are required, but if
52 # the simple dot product in the above algorithm is used then 32 loads are
53 # required. increasing N1 decreases the total number of loads, but only as long
54 # as we have enough temporary registers to keep the matrix blocks and vectors.
55 #
56 #                L                     *    X     =     B
57 #
58 #  [ .                               ]   [ e e ]     [ . . ]
59 #  [ . .                             ]   [ f f ]     [ . . ]
60 #  [ . . .                           ]   [ g g ]     [ . . ]
61 #  [ . . . .                         ]   [ h h ]     [ . . ]
62 #  [ a b c d .                       ]   [ A B ]  =  [ . . ]
63 #  [ a b c d . .                     ]   [ C D ]     [ . . ]
64 #  [ . . . . . . .                   ]   [ . . ]     [ . . ]
65 #  [ . . . . . . . .                 ]   [ . . ]     [ . . ]
66 #  [ . . . . . . . . .               ]   [ . . ]     [ . . ]
67 #
68 # note that L is stored by rows but X and B are stored by columns.
69 # the outer product loops are unrolled for extra speed.
70 #
71 # LDLT FACTORIZATION
72 # ------------------
73 #
74 # the factorization algorithm builds L incrementally by repeatedly solving
75 # the following equation:
76 #
77 #   [ L  0 ] [ D 0 ] [ L' l  ] = [ A  a ]    <-- n rows
78 #   [ l' e ] [ 0 d ] [ 0  e' ]   [ a' b ]    <-- m rows
79 #
80 #   [ L*D*L'   L*D*l         ] = [ A  a ]
81 #   [ l'*D*L'  l'*D*l+e*d*e' ]   [ a' b ]
82 #
83 # L*D*L'=A is an existing solution, and a,b are new rows/columns to add to A.
84 # we compute:
85 #
86 #   L * (Dl) = a
87 #   l = inv(D) * Dl
88 #   e*d*e' = b - l'*Dl       (m*m LDLT factorization)
89 #
90 #
91 # L-transpose solver
92 # ------------------
93 #
94 # the LT (L-transpose) solver uses the same logic as the standard L-solver,
95 # with a few tricks to make it work. to solve L^T*X=B we first remap:
96 #    L to Lhat : Lhat(i,j) = L(n-j,n-i)
97 #    X to Xhat : Xhat(i) = X(n-i)
98 #    B to Bhat : Bhat(i) = B(n-i)
99 # and then solve Lhat*Xhat = Bhat. the current LT solver only supports one
100 # right hand side, but that's okay as it is not used in the factorizer.
101 #
102 #############################################################################
103 #
104 # code generation parameters, set in a parameters file:
105 #    FNAME   : name of source file to generate - a .c file will be made
106 #    TYPE    : 'f' to build factorizer, 's' to build solver, 't' to build the
107 #              transpose solver.
108 #    N1      : block size (size of outer product matrix) (1..9)
109 #    UNROLL1 : solver inner loop unrolling factor (1..)
110 #    UNROLL2 : factorizer inner loop unrolling factor (1..)
111 #    MADD    : if nonzero, generate code for fused multiply-add (0,1)
112 #    FETCH   : how to fetch data in the inner loop:
113 #                0 - load in a batch (the `normal way')
114 #                1 - delay inner loop loads until just before they're needed
115 #
116 #############################################################################
117 #
118 # TODO
119 # ----
120 #
121 # * dFactorLDLT() is not so efficient for matrix sizes < block size, e.g.
122 #   redundant calls, zero loads, adds etc
123 #
124 #############################################################################
125 #
126 # NOTES:
127 #
128 # * on the pentium we can prefetch like this:
129 #     asm ("prefetcht0 %0" : : "m" (*Ai) );
130 #   but it doesn't seem to help much
131
132 require ("BuildUtil");
133
134 # get and check code generation parameters
135 error ("Usage: BuildLDLT <parameters-file>") if $#ARGV != 0;
136 do $ARGV[0];
137
138 if (!defined($FNAME) || !defined($TYPE) || !defined($N1) ||
139     !defined($UNROLL1) || !defined($UNROLL2) || !defined($MADD) ||
140     !defined($FETCH)) {
141   error ("code generation parameters not defined");
142 }
143
144 # check parameters
145 error ("bad TYPE") if $TYPE ne 'f' && $TYPE ne 's' && $TYPE ne 't';
146 error ("bad N1") if $N1 < 1 || $N1 > 9;
147 error ("bad UNROLL1") if $UNROLL1 < 1;
148 error ("bad UNROLL2") if $UNROLL2 < 1;
149 error ("bad MADD") if $MADD != 0 && $MADD != 1;
150 error ("bad FETCH") if $FETCH < 0 && $FETCH > 1;
151
152 #############################################################################
153 # utility
154
155 # functions to handle delayed loading of p and q values.
156 # bit in the the `ploaded' and `qloaded' numbers record what has been loaded,
157 # so we dont load it again.
158
159 sub newLoads
160 {
161   # bits in these numbers say what registers p and q have been loaded so far
162   $ploaded = 0;
163   $qloaded = 0;
164 }
165
166 sub loadedEverything
167 {
168   $ploaded = 0xffffffff;
169   $qloaded = 0xffffffff;
170 }
171
172 sub loadP # (i,loadcmd)
173 {
174   my $i = $_[0];
175   my $loadcmd = $_[1];
176   return if ($ploaded & (1 << $i));
177   output ($loadcmd);
178   $ploaded |= (1 << $i);
179 }
180
181 sub loadQ # (i,loadcmd)
182 {
183   my $i = $_[0];
184   my $loadcmd = $_[1];
185   return if ($qloaded & (1 << $i));
186   output ($loadcmd);
187   $qloaded |= (1 << $i);
188 }
189
190 #############################################################################
191 # make a fast L solve function.
192 # this function has a restriction that the leading dimension of X and B must
193 # be a multiple of the block size.
194
195 sub innerOuterProductLoop # (M,k,nrhs,increment)
196 {
197   my $M=$_[0];
198   my $k=$_[1];
199   my $nrhs=$_[2];
200   my $increment=$_[3];
201   my ($i,$j);
202   newLoads;
203   if ($FETCH==0) {
204     comment ("load p and q values");
205     for ($i=1; $i<=$M; $i++) {
206       if ($TYPE eq 't') {
207         output ("p$i=ell[".ofs2(-($i-1),0,'lskip')."];\n");
208         output ("q$i=ex[".ofs2(-($k),$i-1,'lskip')."];\n") if $i <= $nrhs;
209       }
210       else {
211         output ("p$i=ell[".ofs2($k,$i-1,'lskip')."];\n");
212         output ("q$i=ex[".ofs2($k,$i-1,'lskip')."];\n") if $i <= $nrhs;
213       }
214     }
215     loadedEverything;
216   }
217
218   comment ("compute outer product and add it to the Z matrix");
219   for ($i=1; $i<=$M; $i++) {
220     for ($j=1; $j<=$nrhs; $j++) {
221       if ($TYPE eq 't') {
222         loadP ($i,"p$i=ell[".ofs2(-($i-1),0,'lskip')."];\n");
223         loadQ ($j,"q$j=ex[".ofs2(-($k),$j-1,'lskip')."];\n");
224       }
225       else {
226         loadP ($i,"p$i=ell[".ofs2($k,$i-1,'lskip')."];\n");
227         loadQ ($j,"q$j=ex[".ofs2($k,$j-1,'lskip')."];\n");
228       }
229       my $var = $MADD ? "Z$i$j +=" : "m$i$j =";
230       output ("$var p$i * q$j;\n");
231     }
232   }
233
234   if ($TYPE eq 't') {
235     if ($increment > 0) {
236       output ("ell += lskip1;\n");
237       output ("ex -= $increment;\n");
238     }
239     else {
240       output ("ell += lskip1;\n");
241     }
242   }
243   else {
244     if ($increment > 0) {
245       comment ("advance pointers");
246       output ("ell += $increment;\n");
247       output ("ex += $increment;\n");
248     }
249   }
250
251   if ($MADD==0) {
252     for ($i=1; $i<=$M; $i++) {
253       for ($j=1; $j<=$nrhs; $j++) {
254         output ("Z$i$j += m$i$j;\n");
255       }
256     }
257   }
258 }
259
260
261 sub computeRows # (nrhs,rows)
262 {
263   my $nrhs = $_[0];
264   my $rows = $_[1];
265   my ($i,$j,$k);
266
267   comment ("compute all $rows x $nrhs block of X, from rows i..i+$rows-1");
268
269   comment ("set the Z matrix to 0");
270   for ($i=1; $i<=$rows; $i++) {
271     for ($j=1; $j<=$nrhs; $j++) {
272       output ("Z$i$j=0;\n");
273     }
274   }
275   if ($TYPE eq 't') {
276     output ("ell = L - i;\n");
277   }
278   else {
279     output ("ell = L + i*lskip1;\n");
280   }
281   output ("ex = B;\n");
282
283   comment ("the inner loop that computes outer products and adds them to Z");
284   output ("for (j=i-$UNROLL1; j >= 0; j -= $UNROLL1) {\n");
285   for ($k=0; $k < $UNROLL1; $k++) {
286     innerOuterProductLoop ($rows,$k,$nrhs,($k==$UNROLL1-1) ? $UNROLL1 : 0);
287   }
288
289   comment ("end of inner loop");
290   output ("}\n");
291
292   if ($UNROLL1 > 1) {
293     comment ("compute left-over iterations");
294     output ("j += $UNROLL1;\n");
295     output ("for (; j > 0; j--) {\n");
296     innerOuterProductLoop ($rows,'0',$nrhs,1);
297     output ("}\n");
298   }
299
300   comment ("finish computing the X(i) block");
301
302   for ($j=1; $j<=$nrhs; $j++) {
303     if ($TYPE eq 't') {
304       output ("Z1$j = ex[".ofs1(-($j-1),'lskip')."] - Z1$j;\n");
305       output ("ex[".ofs1(-($j-1),'lskip')."] = Z1$j;\n");
306     }
307     else {
308       output ("Z1$j = ex[".ofs1($j-1,'lskip')."] - Z1$j;\n");
309       output ("ex[".ofs1($j-1,'lskip')."] = Z1$j;\n");
310     }
311   }
312
313   for ($i=2; $i<=$rows; $i++) {
314     for ($j=1; $j<$i; $j++) {
315       if ($TYPE eq 't') {
316         output ("p$j = ell[".ofs2(-($i-1),$j-1,'lskip')."];\n");
317       }
318       else {
319         output ("p$j = ell[".ofs2($j-1,$i-1,'lskip')."];\n");
320       }
321     }
322     for ($j=1; $j<=$nrhs; $j++) {
323       if ($TYPE eq 't') {
324         output ("Z$i$j = ex[".ofs2(-($i-1),$j-1,'lskip')."] - Z$i$j");
325       }
326       else {
327         output ("Z$i$j = ex[".ofs2($i-1,$j-1,'lskip')."] - Z$i$j");
328       }
329       for ($k=1; $k < $i; $k++) {
330         output (" - p$k*Z$k$j");
331       }
332       output (";\n");
333       if ($TYPE eq 't') {
334         output ("ex[".ofs2(-($i-1),$j-1,'lskip')."] = Z$i$j;\n");
335       }
336       else {
337         output ("ex[".ofs2($i-1,$j-1,'lskip')."] = Z$i$j;\n");
338       }
339     }
340   }
341 }
342
343
344 sub makeFastL1Solve # ( number-of-right-hand-sides )
345 {
346   my $nrhs = $_[0];
347   my ($i,$j,$k);
348   my $funcsuffix = ($TYPE eq 'f') ? "_$nrhs" : '';
349   my $staticvoid = ($TYPE eq 'f') ? 'static void' : 'void';
350
351   # function header
352   if ($TYPE eq 't') {
353     output (<<END);
354
355 /* solve L^T * x=b, with b containing 1 right hand side.
356  * L is an n*n lower triangular matrix with ones on the diagonal.
357  * L is stored by rows and its leading dimension is lskip.
358  * b is an n*1 matrix that contains the right hand side.
359  * b is overwritten with x.
360  * this processes blocks of $N1.
361  */
362
363 void dSolveL1T (const dReal *L, dReal *B, int n, int lskip1)
364 {
365 END
366   }
367   else {
368     output (<<END);
369
370 /* solve L*X=B, with B containing $nrhs right hand sides.
371  * L is an n*n lower triangular matrix with ones on the diagonal.
372  * L is stored by rows and its leading dimension is lskip.
373  * B is an n*$nrhs matrix that contains the right hand sides.
374  * B is stored by columns and its leading dimension is also lskip.
375  * B is overwritten with X.
376  * this processes blocks of $N1*$N1.
377  * if this is in the factorizer source file, n must be a multiple of $N1.
378  */
379
380 $staticvoid dSolveL1$funcsuffix (const dReal *L, dReal *B, int n, int lskip1)
381 {
382 END
383   }
384
385   comment ("declare variables - Z matrix, p and q vectors, etc");
386   output ("dReal ");
387   for ($i=1; $i<=$N1; $i++) {
388     for ($j=1; $j<=$nrhs; $j++) {
389       output ("Z$i$j,");                # Z matrix
390       output ("m$i$j,") if ! $MADD;     # temporary vars where multiplies put
391     }
392   }
393   for ($i=1; $i<=$N1; $i++) {
394     output ("p$i,");
395     output ("q$i,") if $i <= $nrhs;
396   }
397   output ("*ex;\nconst dReal *ell;\n");
398   output ("int ");
399   for ($i=2; $i<$N1; $i++) {
400     output ("lskip$i,");
401   }
402   output ("i,j;\n");
403
404   if ($TYPE eq 't') {
405     comment ("special handling for L and B because we're solving L1 *transpose*");
406     output ("L = L + (n-1)*(lskip1+1);\n");
407     output ("B = B + n-1;\n");
408     output ("lskip1 = -lskip1;\n");
409   }
410
411   if ($N1 > 2) {
412     comment ("compute lskip values");
413     for ($i=2; $i<$N1; $i++) {
414       output ("lskip$i = $i*lskip1;\n");
415     }
416   }
417
418   comment ("compute all $N1 x $nrhs blocks of X");
419   if ($TYPE eq 's' or $TYPE eq 't') {
420     output ("for (i=0; i <= n-$N1; i+=$N1) {\n");
421   }
422   else {
423     output ("for (i=0; i < n; i+=$N1) {\n");
424   }
425   computeRows ($nrhs,$N1);
426   comment ("end of outer loop");
427   output ("}\n");
428
429   if ($TYPE eq 's' or $TYPE eq 't') {
430     comment ("compute rows at end that are not a multiple of block size");
431     output ("for (; i < n; i++) {\n");
432     computeRows ($nrhs,1);
433     output ("}\n");
434   }
435
436   output ("}\n");
437 }
438
439 #############################################################################
440 # make a fast L*D*L' factorizer
441
442 # code fragment: this factors an M x M block. if A_or_Z is 0 then it works
443 # on the $A matrix otherwise it works on the Z matrix. in either case it
444 # writes the diagonal entries into the `dee' vector.
445 # it is a simple implementation of the LDLT algorithm, with no tricks.
446
447 sub getA # (i,j,A,A_or_Z)
448 {
449   my $i = $_[0];
450   my $j = $_[1];
451   my $A = $_[2];
452   return $_[3] ? ('Z'.($i+1).($j+1)) : ($A.'['.ofs2($j,$i,'nskip').']');
453 }
454
455 sub miniLDLT # (A,A_or_Z,M)
456 {
457   my ($i,$j,$k);
458   my $A = $_[0];
459   my $AZ = $_[1];
460   my $M = $_[2];
461   comment ("factorize $M x $M block " . ($AZ ? "Z,dee" : "$A,dee"));
462   comment ("factorize row 1");
463   output ("dee[0] = dRecip(".getA(0,0,$A,$AZ).");\n");
464   for ($i=1; $i<$M; $i++) {
465     comment ("factorize row ".($i+1));
466     for ($j=1; $j<$i; $j++) {
467       output (getA($i,$j,$A,$AZ)." -= ");
468       for ($k=0; $k<$j; $k++) {
469         output (" + ") if $k > 0;
470         output (getA($i,$k,$A,$AZ)."*".getA($j,$k,$A,$AZ));
471       }
472       output (";\n");
473     }
474     output ("sum = 0;\n");
475     for ($j=0; $j<$i; $j++) {
476       output ("q1 = ".getA($i,$j,$A,$AZ).";\n");
477       output ("q2 = q1 * dee[$j];\n");
478       output (getA($i,$j,$A,$AZ)." = q2;\n");
479       output ("sum += q1*q2;\n");
480     }
481     output ("dee[$i] = dRecip(".getA($i,$i,$A,$AZ)." - sum);\n");
482   }
483   comment ("done factorizing $M x $M block");
484 }
485
486
487 sub innerScaleAndOuterProductLoop # (M,k)
488 {
489   my $M = $_[0];
490   my $k = $_[1];
491   my ($i,$j);
492   for ($i=1; $i<=$M; $i++) {
493     output ("p$i = ell[".ofs2($k,$i-1,'nskip')."];\n");
494   }
495   output ("dd = dee[$k];\n");
496   for ($i=1; $i<=$M; $i++) {
497     output ("q$i = p$i*dd;\n");
498   }
499   for ($i=1; $i<=$M; $i++) {
500     output ("ell[".ofs2($k,$i-1,'nskip')."] = q$i;\n");
501   }
502   for ($i=1; $i<=$M; $i++) {
503     for ($j=1; $j<=$i; $j++) {
504       my $var = $MADD ? "Z$i$j +=" : "m$i$j =";
505       output ("$var p$i*q$j;\n");
506     }
507   }
508   if ($MADD==0) {
509   for ($i=1; $i<=$M; $i++) {
510     for ($j=1; $j<=$i; $j++) {
511         output ("Z$i$j += m$i$j;\n");
512       }
513     }
514   }
515 }
516
517
518 sub diagRows # (M)
519 {
520   my $M=$_[0];
521   comment ("scale the elements in a $M x i block at A(i,0), and also");
522   comment ("compute Z = the outer product matrix that we'll need.");
523   for ($i=1; $i<=$M; $i++) {
524     for ($j=1; $j<=$i; $j++) {
525       output ("Z$i$j = 0;\n");
526     }
527   }
528   output ("ell = A+i*nskip1;\n");
529   output ("dee = d;\n");
530   output ("for (j=i-$UNROLL2; j >= 0; j -= $UNROLL2) {\n");
531   for ($i=0; $i < $UNROLL2; $i++) {
532     innerScaleAndOuterProductLoop ($M,$i);
533   }
534   output ("ell += $UNROLL2;\n");
535   output ("dee += $UNROLL2;\n");
536   output ("}\n");
537
538   if ($UNROLL2 > 1) {
539     comment ("compute left-over iterations");
540     output ("j += $UNROLL2;\n");
541     output ("for (; j > 0; j--) {\n");
542     innerScaleAndOuterProductLoop ($M,0);
543     output ("ell++;\n");
544     output ("dee++;\n");
545     output ("}\n");
546   }
547 }
548
549
550 sub diagBlock # (M)
551 {
552   my $M = $_[0];
553   comment ("solve for diagonal $M x $M block at A(i,i)");
554   for ($i=1; $i<=$M; $i++) {
555     for ($j=1; $j<=$i; $j++) {
556       output ("Z$i$j = ell[".ofs2($j-1,$i-1,'nskip')."] - Z$i$j;\n");
557     }
558   }
559   output ("dee = d + i;\n");
560   miniLDLT ('',1,$M);
561   for ($i=2; $i<=$M; $i++) {
562     for ($j=1; $j<$i; $j++) {
563       output ("ell[".ofs2($j-1,$i-1,'nskip')."] = Z$i$j;\n");
564     }
565   }
566 }
567
568
569 sub makeFastLDLT
570 {
571   my ($i,$j,$k);
572
573   # function header
574   output (<<END);
575
576
577 void dFactorLDLT (dReal *A, dReal *d, int n, int nskip1)
578 {
579 END
580   output ("int i,j");
581   for ($i=2; $i<$N1; $i++) {
582     output (",nskip$i");
583   }
584   output (";\n");
585   output ("dReal sum,*ell,*dee,dd,p1,p2");
586   for ($i=3; $i<=$N1; $i++) {
587     output (",p$i");
588   }
589   for ($i=1; $i<=$N1; $i++) {
590     output (",q$i");
591   }
592   for ($i=1; $i<=$N1; $i++) {
593     for ($j=1; $j<=$i; $j++) {
594       output (",Z$i$j");
595       output (",m$i$j") if ! $MADD;     # temporary vars where multiplies put
596     }
597   }
598   output (";\n");
599   output ("if (n < 1) return;\n");
600   # output ("nskip1 = dPAD(n);\n");   ... not any more
601   for ($i=2; $i<$N1; $i++) {
602     output ("nskip$i = $i*nskip1;\n");
603   }
604
605   output ("\nfor (i=0; i<=n-$N1; i += $N1) {\n");
606   comment ("solve L*(D*l)=a, l is scaled elements in $N1 x i block at A(i,0)");
607   output ("dSolveL1_$N1 (A,A+i*nskip1,i,nskip1);\n");
608
609   diagRows ($N1);
610   diagBlock ($N1);
611   output ("}\n");
612
613   comment ("compute the (less than $N1) rows at the bottom");
614   output ("switch (n-i) {\n");
615   output ("case 0:\n");
616   output ("break;\n\n");
617
618   for ($i=1; $i<$N1; $i++) {
619     output ("case $i:\n");
620     output ("dSolveL1_$i (A,A+i*nskip1,i,nskip1);\n");
621     diagRows ($i);
622     diagBlock ($i);
623     output ("break;\n\n");
624   }
625
626   output ("default: *((char*)0)=0;  /* this should never happen! */\n");
627   output ("}\n");
628
629   output ("}\n");
630 }
631
632 #############################################################################
633 # write source code
634
635 open (FOUT,">$FNAME.c") or die "can't open $FNAME.c for writing";
636
637 # file and function header
638 output (<<END);
639 /* generated code, do not edit. */
640
641 #include "ode/matrix.h"
642 END
643
644 if ($TYPE eq 'f') {
645   for ($i=1; $i <= $N1; $i++) {
646     makeFastL1Solve ($i);
647   }
648   makeFastLDLT;
649 }
650 else {
651   makeFastL1Solve (1);
652   makeRealFastL1Solve;
653 }
654 close FOUT;