gaussian_elimination_export.cpp
Go to the documentation of this file.
1 /*
2  * This file is part of ACADO Toolkit.
3  *
4  * ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
5  * Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
6  * Milan Vukov, Rien Quirynen, KU Leuven.
7  * Developed within the Optimization in Engineering Center (OPTEC)
8  * under supervision of Moritz Diehl. All rights reserved.
9  *
10  * ACADO Toolkit is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 3 of the License, or (at your option) any later version.
14  *
15  * ACADO Toolkit is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  * Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with ACADO Toolkit; if not, write to the Free Software
22  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23  *
24  */
25 
26 
27 
34 #include <acado/code_generation/linear_solvers/gaussian_elimination_export.hpp>
35 
36 using namespace std;
37 
38 BEGIN_NAMESPACE_ACADO
39 
40 //
41 // PUBLIC MEMBER FUNCTIONS:
42 //
43 
44 ExportGaussElim::ExportGaussElim( UserInteraction* _userInteraction,
45  const std::string& _commonHeaderName
46  ) : ExportLinearSolver( _userInteraction,_commonHeaderName )
47 {
48 }
49 
50 ExportGaussElim::~ExportGaussElim( )
51 {}
52 
53 returnValue ExportGaussElim::getDataDeclarations( ExportStatementBlock& declarations,
54  ExportStruct dataStruct
55  ) const
56 {
57  declarations.addDeclaration( rk_swap,dataStruct ); // needed for the row swaps
58  if( REUSE ) {
59  declarations.addDeclaration( rk_bPerm,dataStruct ); // reordered right-hand side
60  if( TRANSPOSE ) {
61  declarations.addDeclaration( rk_bPerm_trans,dataStruct );
62  }
63  }
64 
65  return SUCCESSFUL_RETURN;
66 }
67 
68 
69 returnValue ExportGaussElim::getFunctionDeclarations( ExportStatementBlock& declarations
70  ) const
71 {
72  declarations.addDeclaration( solve );
73  declarations.addDeclaration( solveTriangular );
74  if( REUSE ) {
75  declarations.addDeclaration( solveReuse );
76  }
77 
78  return SUCCESSFUL_RETURN;
79 }
80 
81 
82 returnValue ExportGaussElim::getCode( ExportStatementBlock& code
83  )
84 {
85 
86  if (nRightHandSides > 0) {
87  if( !REUSE ) return ACADOERROR(RET_INVALID_OPTION);
88 
89  setupFactorization( solve, rk_swap, determinant, string("fabs") );
90  code.addFunction( solve );
91 
92  setupSolveReuseComplete( solveReuse, rk_bPerm );
93  code.addFunction( solveReuse );
94 
95  if( TRANSPOSE ) {
96  setupSolveReuseTranspose( solveReuseTranspose, rk_bPerm_trans );
97  code.addFunction( solveReuseTranspose );
98  }
99  }
100  else {
101  setupSolveUpperTriangular( solveTriangular );
102  code.addFunction( solveTriangular );
103 
104  setupSolve( solve, solveTriangular, rk_swap, determinant, string("fabs") );
105  code.addFunction( solve );
106 
107  if( REUSE ) { // Also export the extra function which reuses the factorization of the matrix A
108 
109  setupSolveReuse( solveReuse, solveTriangular, rk_bPerm );
110  code.addFunction( solveReuse );
111  }
112 // if( TRANSPOSE ) return ACADOERROR( RET_NOT_YET_IMPLEMENTED );
113  if( TRANSPOSE ) {
114  setupSolveReuseTranspose( solveReuseTranspose, rk_bPerm_trans );
115  code.addFunction( solveReuseTranspose );
116  }
117  }
118 
119  return SUCCESSFUL_RETURN;
120 }
121 
122 
123 returnValue ExportGaussElim::setupSolveUpperTriangular( ExportFunction& _solveTriangular ) {
124 
125  if (nRightHandSides > 0)
126  return ACADOERROR(RET_INVALID_OPTION);
127 
128  uint run1, run2;
129  // Solve the upper triangular system of equations:
130  for( run1 = dim; run1 > 0; run1--) {
131  for( run2 = dim-1; run2 > (run1-1); run2--) {
132  _solveTriangular.addStatement( b.getRow( (run1-1) ) -= A.getSubMatrix( (run1-1),(run1-1)+1,run2,run2+1 ) * b.getRow( run2 ) );
133  }
134  _solveTriangular << "b[" << toString( (run1-1) ) << "] = b[" << toString( (run1-1) ) << "]/A[" << toString( (run1-1)*dim+(run1-1) ) << "];\n";
135  }
136 
137  return SUCCESSFUL_RETURN;
138 }
139 
140 
141 returnValue ExportGaussElim::setupSolve( ExportFunction& _solve, ExportFunction& _solveTriangular, ExportVariable& _swap, ExportVariable& _determinant, const string& absF ) {
142 
143  uint run1, run2, run3;
144  ExportIndex i( "i" );
145  _solve.addIndex( i );
146  ExportIndex j( "j" );
147  ExportIndex k( "k" );
148  ExportVariable indexMax( "indexMax", 1, 1, INT, ACADO_LOCAL, true );
149  ExportVariable intSwap( "intSwap", 1, 1, INT, ACADO_LOCAL, true );
150  ExportVariable valueMax( "valueMax", 1, 1, REAL, ACADO_LOCAL, true );
151  ExportVariable temp( "temp", 1, 1, REAL, ACADO_LOCAL, true );
152  if( !UNROLLING ) {
153  _solve.addIndex( j );
154  _solve.addIndex( k );
155  _solve.addDeclaration( indexMax );
156  if( REUSE ) _solve.addDeclaration( intSwap );
157  _solve.addDeclaration( valueMax );
158  _solve.addDeclaration( temp );
159  }
160 
161  if (nRightHandSides > 0)
162  return ACADOERROR(RET_INVALID_OPTION);
163 
164  // initialise rk_perm (the permutation vector)
165  if( REUSE ) {
166  ExportForLoop loop1( i,0,dim );
167  loop1 << rk_perm.get( 0,i ) << " = " << i.getName() << ";\n";
168  _solve.addStatement( loop1 );
169  }
170 
171  _solve.addStatement( _determinant == 1 );
172  if( UNROLLING || dim <= 5 ) {
173  // Start the factorization:
174  for( run1 = 0; run1 < (dim-1); run1++ ) {
175  // Search for pivot in column run1:
176  for( run2 = run1+1; run2 < dim; run2++ ) {
177  // add the test (if or else if):
178  stringstream test;
179  if( run2 == (run1+1) ) {
180  test << "if(";
181  } else {
182  test << "else if(";
183  }
184  test << absF << "(A[" << toString( run2*dim+run1 ) << "]) > " << absF << "(A[" << toString( run1*dim+run1 ) << "])";
185  for( run3 = run1+1; run3 < dim; run3++ ) {
186  if( run3 != run2) {
187  test << " && " << absF << "(A[" << toString( run2*dim+run1 ) << "]) > " << absF << "(A[" << toString( run3*dim+run1 ) << "])";
188  }
189  }
190  test << ") {\n";
191  _solve.addStatement( test.str() );
192 
193  // do the row swaps:
194  // for A:
195  for( run3 = 0; run3 < dim; run3++ ) {
196  _solve.addStatement( _swap == A.getSubMatrix( run1,run1+1,run3,run3+1 ) );
197  _solve.addStatement( A.getSubMatrix( run1,run1+1,run3,run3+1 ) == A.getSubMatrix( run2,run2+1,run3,run3+1 ) );
198  _solve.addStatement( A.getSubMatrix( run2,run2+1,run3,run3+1 ) == _swap );
199  }
200  // for b:
201  _solve.addStatement( _swap == b.getRow( run1 ) );
202  _solve.addStatement( b.getRow( run1 ) == b.getRow( run2 ) );
203  _solve.addStatement( b.getRow( run2 ) == _swap );
204 
205  if( REUSE ) { // rk_perm also needs to be updated if it needs to be possible to reuse the factorization
206  _solve.addStatement( intSwap == rk_perm.getCol( run1 ) );
207  _solve.addStatement( rk_perm.getCol( run1 ) == rk_perm.getCol( run2 ) );
208  _solve.addStatement( rk_perm.getCol( run2 ) == intSwap );
209  }
210 
211  _solve.addStatement( "}\n" );
212  }
213  // potentially needed row swaps are done
214  _solve.addLinebreak();
215  // update of the next rows:
216  for( run2 = run1+1; run2 < dim; run2++ ) {
217  _solve << "A[" << toString( run2*dim+run1 ) << "] = -A[" << toString( run2*dim+run1 ) << "]/A[" << toString( run1*dim+run1 ) << "];\n";
218  _solve.addStatement( A.getSubMatrix( run2,run2+1,run1+1,dim ) += A.getSubMatrix( run2,run2+1,run1,run1+1 ) * A.getSubMatrix( run1,run1+1,run1+1,dim ) );
219  _solve.addStatement( b.getRow( run2 ) += A.getSubMatrix( run2,run2+1,run1,run1+1 ) * b.getRow( run1 ) );
220  _solve.addLinebreak();
221  }
222  _solve.addStatement( _determinant == _determinant*A.getSubMatrix(run1,run1+1,run1,run1+1) );
223  _solve.addLinebreak();
224  }
225  _solve.addStatement( _determinant == _determinant*A.getSubMatrix(dim-1,dim,dim-1,dim) );
226  _solve.addLinebreak();
227  }
228  else { // without UNROLLING:
229  _solve << "for( i=0; i < (" << toString( dim-1 ) << "); i++ ) {\n";
230  _solve << " indexMax = i;\n";
231  _solve << " valueMax = " << absF << "(A[i*" << toString( dim ) << "+i]);\n";
232  _solve << " for( j=(i+1); j < " << toString( dim ) << "; j++ ) {\n";
233  _solve << " temp = " << absF << "(A[j*" << toString( dim ) << "+i]);\n";
234  _solve << " if( temp > valueMax ) {\n";
235  _solve << " indexMax = j;\n";
236  _solve << " valueMax = temp;\n";
237  _solve << " }\n";
238  _solve << " }\n";
239  _solve << " if( indexMax > i ) {\n";
240  ExportForLoop loop2( k,0,dim );
241  loop2 << " " << _swap.getFullName() << " = A[i*" << toString( dim ) << "+" << k.getName() << "];\n";
242  loop2 << " A[i*" << toString( dim ) << "+" << k.getName() << "] = A[indexMax*" << toString( dim ) << "+" << k.getName() << "];\n";
243  loop2 << " A[indexMax*" << toString( dim ) << "+" << k.getName() << "] = " << _swap.getFullName() << ";\n";
244  _solve.addStatement( loop2 );
245  _solve << " " << _swap.getFullName() << " = b[i];\n";
246  _solve << " b[i] = b[indexMax];\n";
247  _solve << " b[indexMax] = " << _swap.getFullName() << ";\n";
248  if( REUSE ) {
249  _solve << " " << intSwap.getFullName() << " = " << rk_perm.getFullName() << "[i];\n";
250  _solve << " " << rk_perm.getFullName() << "[i] = " << rk_perm.getFullName() << "[indexMax];\n";
251  _solve << " " << rk_perm.getFullName() << "[indexMax] = " << intSwap.getFullName() << ";\n";
252  }
253  _solve << " }\n";
254  _solve << " " << _determinant.getFullName() << " *= A[i*" << toString( dim ) << "+i];\n";
255  _solve << " for( j=i+1; j < " << toString( dim ) << "; j++ ) {\n";
256  _solve << " A[j*" << toString( dim ) << "+i] = -A[j*" << toString( dim ) << "+i]/A[i*" << toString( dim ) << "+i];\n";
257  _solve << " for( k=i+1; k < " << toString( dim ) << "; k++ ) {\n";
258  _solve << " A[j*" << toString( dim ) << "+k] += A[j*" << toString( dim ) << "+i] * A[i*" << toString( dim ) << "+k];\n";
259  _solve << " }\n";
260  _solve << " b[j] += A[j*" << toString( dim ) << "+i] * b[i];\n";
261  _solve << " }\n";
262  _solve << "}\n";
263  _solve << _determinant.getFullName() << " *= A[" << toString( (dim-1)*dim+(dim-1) ) << "];\n";
264  }
265  _solve << _determinant.getFullName() << " = " << absF << "(" << _determinant.getFullName() << ");\n";
266 
267  _solve.addFunctionCall( _solveTriangular, A, b );
268 
269  return SUCCESSFUL_RETURN;
270 }
271 
272 
273 returnValue ExportGaussElim::setupFactorization( ExportFunction& _solve, ExportVariable& _swap, ExportVariable& _determinant, const string& absF ) {
274 
275  uint run1, run2, run3;
276  ExportIndex i( "i" );
277  _solve.addIndex( i );
278  ExportIndex j( "j" );
279  ExportIndex k( "k" );
280  ExportVariable indexMax( "indexMax", 1, 1, INT, ACADO_LOCAL, true );
281  ExportVariable intSwap( "intSwap", 1, 1, INT, ACADO_LOCAL, true );
282  ExportVariable valueMax( "valueMax", 1, 1, REAL, ACADO_LOCAL, true );
283  ExportVariable temp( "temp", 1, 1, REAL, ACADO_LOCAL, true );
284  if( !UNROLLING ) {
285  _solve.addIndex( j );
286  _solve.addIndex( k );
287  _solve.addDeclaration( indexMax );
288  if( REUSE ) _solve.addDeclaration( intSwap );
289  _solve.addDeclaration( valueMax );
290  _solve.addDeclaration( temp );
291  }
292 
293  // initialise rk_perm (the permutation vector)
294  if( REUSE ) {
295  ExportForLoop loop1( i,0,dim );
296  loop1 << rk_perm.get( 0,i ) << " = " << i.getName() << ";\n";
297  _solve.addStatement( loop1 );
298  }
299 
300  _solve.addStatement( _determinant == 1 );
301  if( UNROLLING || dim <= 5 ) {
302  // Start the factorization:
303  for( run1 = 0; run1 < (dim-1); run1++ ) {
304  // Search for pivot in column run1:
305  for( run2 = run1+1; run2 < dim; run2++ ) {
306  // add the test (if or else if):
307  stringstream test;
308  if( run2 == (run1+1) ) {
309  test << "if(";
310  } else {
311  test << "else if(";
312  }
313  test << absF << "(A[" << toString( run2*dim+run1 ) << "]) > " << absF << "(A[" << toString( run1*dim+run1 ) << "])";
314  for( run3 = run1+1; run3 < dim; run3++ ) {
315  if( run3 != run2) {
316  test << " && " << absF << "(A[" << toString( run2*dim+run1 ) << "]) > " << absF << "(A[" << toString( run3*dim+run1 ) << "])";
317  }
318  }
319  test << ") {\n";
320  _solve.addStatement( test.str() );
321 
322  // do the row swaps:
323  // for A:
324  for( run3 = 0; run3 < dim; run3++ ) {
325  _solve.addStatement( _swap == A.getSubMatrix( run1,run1+1,run3,run3+1 ) );
326  _solve.addStatement( A.getSubMatrix( run1,run1+1,run3,run3+1 ) == A.getSubMatrix( run2,run2+1,run3,run3+1 ) );
327  _solve.addStatement( A.getSubMatrix( run2,run2+1,run3,run3+1 ) == _swap );
328  }
329 
330  if( REUSE ) { // rk_perm also needs to be updated if it needs to be possible to reuse the factorization
331  _solve.addStatement( intSwap == rk_perm.getCol( run1 ) );
332  _solve.addStatement( rk_perm.getCol( run1 ) == rk_perm.getCol( run2 ) );
333  _solve.addStatement( rk_perm.getCol( run2 ) == intSwap );
334  }
335 
336  _solve.addStatement( "}\n" );
337  }
338  // potentially needed row swaps are done
339  _solve.addLinebreak();
340  // update of the next rows:
341  for( run2 = run1+1; run2 < dim; run2++ ) {
342  _solve << "A[" << toString( run2*dim+run1 ) << "] = -A[" << toString( run2*dim+run1 ) << "]/A[" << toString( run1*dim+run1 ) << "];\n";
343  _solve.addStatement( A.getSubMatrix( run2,run2+1,run1+1,dim ) += A.getSubMatrix( run2,run2+1,run1,run1+1 ) * A.getSubMatrix( run1,run1+1,run1+1,dim ) );
344  _solve.addLinebreak();
345  }
346  _solve.addStatement( _determinant == _determinant*A.getSubMatrix(run1,run1+1,run1,run1+1) );
347  _solve.addLinebreak();
348  }
349  _solve.addStatement( _determinant == _determinant*A.getSubMatrix(dim-1,dim,dim-1,dim) );
350  _solve.addLinebreak();
351  }
352  else { // without UNROLLING:
353  _solve << "for( i=0; i < (" << toString( dim-1 ) << "); i++ ) {\n";
354  _solve << " indexMax = i;\n";
355  _solve << " valueMax = " << absF << "(A[i*" << toString( dim ) << "+i]);\n";
356  _solve << " for( j=(i+1); j < " << toString( dim ) << "; j++ ) {\n";
357  _solve << " temp = " << absF << "(A[j*" << toString( dim ) << "+i]);\n";
358  _solve << " if( temp > valueMax ) {\n";
359  _solve << " indexMax = j;\n";
360  _solve << " valueMax = temp;\n";
361  _solve << " }\n";
362  _solve << " }\n";
363  _solve << " if( indexMax > i ) {\n";
364  ExportForLoop loop2( k,0,dim );
365  loop2 << " " << _swap.getFullName() << " = A[i*" << toString( dim ) << "+" << k.getName() << "];\n";
366  loop2 << " A[i*" << toString( dim ) << "+" << k.getName() << "] = A[indexMax*" << toString( dim ) << "+" << k.getName() << "];\n";
367  loop2 << " A[indexMax*" << toString( dim ) << "+" << k.getName() << "] = " << _swap.getFullName() << ";\n";
368  _solve.addStatement( loop2 );
369  if( REUSE ) {
370  _solve << " " << intSwap.getFullName() << " = " << rk_perm.getFullName() << "[i];\n";
371  _solve << " " << rk_perm.getFullName() << "[i] = " << rk_perm.getFullName() << "[indexMax];\n";
372  _solve << " " << rk_perm.getFullName() << "[indexMax] = " << intSwap.getFullName() << ";\n";
373  }
374  _solve << " }\n";
375  _solve << " " << _determinant.getFullName() << " *= A[i*" << toString( dim ) << "+i];\n";
376  _solve << " for( j=i+1; j < " << toString( dim ) << "; j++ ) {\n";
377  _solve << " A[j*" << toString( dim ) << "+i] = -A[j*" << toString( dim ) << "+i]/A[i*" << toString( dim ) << "+i];\n";
378  _solve << " for( k=i+1; k < " << toString( dim ) << "; k++ ) {\n";
379  _solve << " A[j*" << toString( dim ) << "+k] += A[j*" << toString( dim ) << "+i] * A[i*" << toString( dim ) << "+k];\n";
380  _solve << " }\n";
381  _solve << " }\n";
382  _solve << "}\n";
383  _solve << _determinant.getFullName() << " *= A[" << toString( (dim-1)*dim+(dim-1) ) << "];\n";
384  }
385  _solve << _determinant.getFullName() << " = " << absF << "(" << _determinant.getFullName() << ");\n";
386 
387  return SUCCESSFUL_RETURN;
388 }
389 
390 
391 returnValue ExportGaussElim::setupSolveReuse( ExportFunction& _solveReuse, ExportFunction& _solveTriangular, ExportVariable& _bPerm ) {
392 
393  uint run1, run2;
394 
395  if (nRightHandSides > 0)
396  return ACADOERROR(RET_INVALID_OPTION);
397 
398  for( run1 = 0; run1 < dim; run1++ ) {
399  _solveReuse << _bPerm.get( run1,0 ) << " = b[" << rk_perm.getFullName() << "[" << toString( run1 ) << "]];\n";
400  }
401 
402  for( run2 = 1; run2 < dim; run2++ ) { // row run2
403  for( run1 = 0; run1 < run2; run1++ ) { // column run1
404  _solveReuse << _bPerm.get( run2,0 ) << " += A[" << toString( run2*dim+run1 ) << "]*" << _bPerm.getFullName() << "[" << toString( run1 ) << "];\n";
405  }
406  _solveReuse.addLinebreak();
407  }
408  _solveReuse.addLinebreak();
409 
410  _solveReuse.addFunctionCall( _solveTriangular, A, _bPerm );
411  _solveReuse.addStatement( b == _bPerm );
412 
413  return SUCCESSFUL_RETURN;
414 }
415 
416 
417 returnValue ExportGaussElim::setupSolveReuseComplete( ExportFunction& _solveReuse, ExportVariable& _bPerm ) {
418 
419  ExportIndex run1( "i" );
420  ExportIndex run2( "j" );
421  ExportIndex tmp_index1( "index1" );
422  ExportIndex tmp_index2( "index2" );
423  ExportVariable tmp( "tmp_var", 1, 1, _bPerm.getType(), ACADO_LOCAL, true );
424  _solveReuse.addIndex( run1 );
425  _solveReuse.addIndex( run2 );
426  _solveReuse.addIndex( tmp_index1 );
427  _solveReuse.addIndex( tmp_index2 );
428  _solveReuse.addDeclaration(tmp);
429  uint run3;
430 
431  if (nRightHandSides <= 0)
432  return ACADOERROR(RET_INVALID_OPTION);
433 
434  ExportForLoop loop1( run1, 0, dim );
435  loop1 << run2.getName() << " = " << rk_perm.getFullName() << "[" << run1.getName() << "]*" << toString(nRightHandSides) << ";\n";
436  for( run3 = 0; run3 < nRightHandSides; run3++ ) {
437  loop1 << _bPerm.get( run1,run3 ) << " = b[" << run2.getName() << "+" << toString(run3) << "];\n";
438  }
439  _solveReuse.addStatement( loop1 );
440 
441  ExportForLoop loop2( run2, 1, dim ); // row run2
442  loop2.addStatement( tmp_index1 == run2*nRightHandSides );
443  ExportForLoop loop3( run1, 0, run2 ); // column run1
444  loop3.addStatement( tmp_index2 == run1*nRightHandSides );
445  loop3.addStatement( tmp == A.getElement(run2,run1) );
446  for( run3 = 0; run3 < nRightHandSides; run3++ ) {
447 // loop3.addStatement( _bPerm.getElement( run2,run3 ) += tmp * _bPerm.getElement( run1,run3 ) );
448  loop3 << _bPerm.getFullName() << "[" << tmp_index1.getName() << "+" << toString(run3) << "] += " << tmp.getName() << "*" << _bPerm.getFullName() << "[" << tmp_index2.getName() << "+" << toString(run3) << "];\n";
449  }
450  loop2.addStatement( loop3 );
451  _solveReuse.addStatement( loop2 );
452 
453 
454  // Solve the upper triangular system of equations:
455  ExportForLoop loop4( run1, dim-1, -1, -1 );
456  loop4.addStatement( tmp_index1 == run1*nRightHandSides );
457  ExportForLoop loop5( run2, dim-1, run1, -1 );
458  loop5.addStatement( tmp_index2 == run2*nRightHandSides );
459  loop5.addStatement( tmp == A.getElement( run1,run2 ) );
460  for( run3 = 0; run3 < nRightHandSides; run3++ ) {
461 // loop5.addStatement( _bPerm.getElement( run1,run3 ) -= tmp * _bPerm.getElement( run2,run3 ) );
462  loop5 << _bPerm.getFullName() << "[" << tmp_index1.getName() << "+" << toString(run3) << "] -= " << tmp.getName() << "*" << _bPerm.getFullName() << "[" << tmp_index2.getName() << "+" << toString(run3) << "];\n";
463  }
464  loop4.addStatement( loop5 );
465  loop4 << tmp.getName() << " = 1.0/A[" << run1.getName() << "*" << toString(dim+1) << "];\n";
466  for( run3 = 0; run3 < nRightHandSides; run3++ ) {
467 // loop4 << _bPerm.get( run1,run3 ) << " = " << _bPerm.get( run1,run3 ) << "*" << tmp.getName() << ";\n";
468  loop4 << _bPerm.getFullName() << "[" << tmp_index1.getName() << "+" << toString(run3) << "] = " << tmp.getName() << "*" << _bPerm.getFullName() << "[" << tmp_index1.getName() << "+" << toString(run3) << "];\n";
469  }
470  _solveReuse.addStatement( loop4 );
471 
472  _solveReuse.addStatement( b == _bPerm );
473 
474  return SUCCESSFUL_RETURN;
475 }
476 
477 
478 returnValue ExportGaussElim::setupSolveReuseTranspose( ExportFunction& _solveReuse, ExportVariable& _bPerm ) {
479 
480  ExportIndex run1( "i" );
481  ExportIndex run2( "j" );
482  ExportVariable tmp( "tmp_var", 1, 1, _bPerm.getType(), ACADO_LOCAL, true );
483  _solveReuse.addIndex( run1 );
484  _solveReuse.addIndex( run2 );
485  _solveReuse.addDeclaration(tmp);
486 
487  _solveReuse.addStatement( _bPerm == b_trans );
488 
489  ExportForLoop loop2( run2, 0, dim ); // row run2
490  ExportForLoop loop3( run1, 0, run2 ); // column run1
491  loop3.addStatement( _bPerm.getRow(run2) -= A.getElement(run1,run2)*_bPerm.getRow(run1) );
492  loop2.addStatement( loop3 );
493  loop2 << tmp.getName() << " = 1.0/A[" << run2.getName() << "*" << toString(dim+1) << "];\n";
494  loop2.addStatement( _bPerm.getRow(run2) == _bPerm.getRow(run2)*tmp );
495  _solveReuse.addStatement( loop2 );
496 
497 
498  // Solve the upper triangular system of equations:
499  ExportForLoop loop4( run1, dim-1, -1, -1 );
500  ExportForLoop loop5( run2, dim-1, run1, -1 );
501  loop5.addStatement( _bPerm.getRow(run1) += A.getElement(run2,run1)*_bPerm.getRow(run2) );
502  loop4.addStatement( loop5 );
503  _solveReuse.addStatement( loop4 );
504 
505 
506  // The permutation now happens HERE!
507  ExportForLoop loop1( run1, 0, dim );
508  loop1 << run2.getName() << " = " << rk_perm.getFullName() << "[" << run1.getName() << "];\n";
509  loop1.addStatement( b_trans.getRow(run2) == _bPerm.getRow(run1) );
510  _solveReuse.addStatement( loop1 );
511 
512  return SUCCESSFUL_RETURN;
513 }
514 
515 
516 returnValue ExportGaussElim::appendVariableNames( stringstream& string ) {
517 
518  string << ", " << rk_swap.getFullName();
519  if( REUSE ) {
520 // string << ", " << rk_perm.getFullName().getName();
521  string << ", " << rk_bPerm.getFullName();
522  }
523 
524  return SUCCESSFUL_RETURN;
525 }
526 
527 
528 returnValue ExportGaussElim::setup( )
529 {
530  // Other cases are not implemented...
531  ASSERT_RETURN(nCols == nRows);
532 
533  int useOMP;
534  get(CG_USE_OPENMP, useOMP);
535  ExportStruct structWspace;
536  structWspace = useOMP ? ACADO_LOCAL : ACADO_WORKSPACE;
537 
538  rk_swap = ExportVariable( std::string( "rk_" ) + identifier + "swap", 1, 1, REAL, structWspace, true );
539  A = ExportVariable( "A", dim, dim, REAL );
540  if (nRightHandSides > 0) {
541  b = ExportVariable( "b", dim, nRightHandSides, REAL );
542  rk_bPerm = ExportVariable( std::string( "rk_" ) + identifier + "bPerm", dim, nRightHandSides, REAL, structWspace );
543  solve = ExportFunction( getNameSolveFunction(), A, rk_perm );
544  }
545  else {
546  b = ExportVariable( "b", dim, 1, REAL );
547  rk_bPerm = ExportVariable( std::string( "rk_" ) + identifier + "bPerm", dim, 1, REAL, structWspace );
548  solve = ExportFunction( getNameSolveFunction(), A, b, rk_perm );
549  }
550  rk_perm = ExportVariable( "rk_perm", 1, dim, INT );
551  solve.setReturnValue( determinant, false );
552  solve.addLinebreak( ); // FIX: TO MAKE SURE IT GETS EXPORTED
553  if (nRightHandSides <= 0) {
554  solveTriangular = ExportFunction( std::string( "solve_" ) + identifier + "triangular", A, b );
555  solveTriangular.addLinebreak( ); // FIX: TO MAKE SURE IT GETS EXPORTED
556  }
557 
558  if( REUSE ) {
559  solveReuse = ExportFunction( getNameSolveReuseFunction(), A, b, rk_perm );
560  solveReuse.addLinebreak( ); // FIX: TO MAKE SURE IT GETS EXPORTED
561  if( TRANSPOSE ) {
562  b_trans = ExportVariable( "b", dim, 1, REAL );
563  rk_bPerm_trans = ExportVariable( std::string( "rk_" ) + identifier + "bPerm_trans", dim, 1, REAL, structWspace );
564  solveReuseTranspose = ExportFunction( getNameSolveTransposeReuseFunction(), A, b, rk_perm );
565  solveReuseTranspose.addLinebreak( ); // FIX: TO MAKE SURE IT GETS EXPORTED
566  }
567  }
568 
569  int unrollOpt;
570  userInteraction->get( UNROLL_LINEAR_SOLVER, unrollOpt );
571  UNROLLING = (bool) unrollOpt;
572 
573  return SUCCESSFUL_RETURN;
574 }
575 
576 
577 ExportVariable ExportGaussElim::getGlobalExportVariable( const uint factor ) const {
578 
579  int useOMP;
580  get(CG_USE_OPENMP, useOMP);
581  ExportStruct structWspace;
582  structWspace = useOMP ? ACADO_LOCAL : ACADO_WORKSPACE;
583 
584  return ExportVariable( std::string( "rk_" ) + identifier + "perm", factor, dim, INT, structWspace );
585 }
586 
587 
588 //
589 // PROTECTED MEMBER FUNCTIONS:
590 //
591 
592 
593 
594 CLOSE_NAMESPACE_ACADO
595 
596 // end of file.
ExportVariable::getRow
ExportVariable getRow(const ExportIndex &idx) const
Definition: export_variable.cpp:395
ExportLinearSolver::getNameSolveTransposeReuseFunction
const std::string getNameSolveTransposeReuseFunction()
Definition: linear_solver_export.cpp:187
ExportLinearSolver::solveReuseTranspose
ExportFunction solveReuseTranspose
Definition: linear_solver_export.hpp:296
SUCCESSFUL_RETURN
Definition: acado_types.hpp:836
ASSERT_RETURN
#define ASSERT_RETURN(x)
Definition: acado_debugging.hpp:52
ACADO_WORKSPACE
Definition: acado_types.hpp:805
AlgorithmicBase::userInteraction
UserInteraction * userInteraction
Definition: algorithmic_base.hpp:597
ExportGaussElim::setup
virtual returnValue setup()
Definition: gaussian_elimination_export.cpp:528
Options::get
returnValue get(OptionsName name, int &value) const
Definition: options.cpp:69
gaussian_elimination_export.hpp
RET_INVALID_OPTION
Definition: acado_types.hpp:889
returnValue
Allows to pass back messages to the calling function.
Definition: acado_types.hpp:1139
uint
BEGIN_NAMESPACE_ACADO typedef unsigned int uint
Definition: acado_types.hpp:42
ExportGaussElim::setupSolveReuseTranspose
virtual returnValue setupSolveReuseTranspose(ExportFunction &_solveReuse, ExportVariable &_bPerm)
Definition: gaussian_elimination_export.cpp:478
ExportForLoop
Allows to export code of a for-loop.
Definition: export_for_loop.hpp:56
std::toString
string toString(T const &value)
Definition: acado_io_utils.hpp:208
ExportGaussElim::setupSolveReuseComplete
virtual returnValue setupSolveReuseComplete(ExportFunction &_solveReuse, ExportVariable &_bPerm)
Definition: gaussian_elimination_export.cpp:417
ExportVariable::getElement
ExportVariable getElement(const ExportIndex &rowIdx, const ExportIndex &colIdx) const
Definition: export_variable.cpp:434
CLOSE_NAMESPACE_ACADO
#define CLOSE_NAMESPACE_ACADO
Definition: acado_namespace_macros.hpp:58
ExportVariable::getSubMatrix
ExportVariable getSubMatrix(const ExportIndex &rowIdx1, const ExportIndex &rowIdx2, const ExportIndex &colIdx1, const ExportIndex &colIdx2) const
Definition: export_variable.cpp:425
ExportGaussElim::getCode
virtual returnValue getCode(ExportStatementBlock &code)
Definition: gaussian_elimination_export.cpp:82
ExportIndex
Defines a scalar-valued index variable to be used for exporting code.
Definition: export_index.hpp:54
ExportLinearSolver::solveTriangular
ExportFunction solveTriangular
Definition: linear_solver_export.hpp:294
ExportGaussElim::appendVariableNames
returnValue appendVariableNames(std::stringstream &string)
Definition: gaussian_elimination_export.cpp:516
ExportStruct
ExportStruct
Definition: acado_types.hpp:802
INT
Definition: acado_types.hpp:795
ExportLinearSolver::getNameSolveReuseFunction
const std::string getNameSolveReuseFunction()
Definition: linear_solver_export.cpp:181
ExportVariable::get
const std::string get(const ExportIndex &rowIdx, const ExportIndex &colIdx) const
Definition: export_variable.cpp:227
UNROLL_LINEAR_SOLVER
Definition: acado_types.hpp:349
UserInteraction
Encapsulates all user interaction for setting options, logging data and plotting results.
Definition: user_interaction.hpp:57
ExportFunction
Allows to export code of an arbitrary function.
Definition: export_function.hpp:58
ExportStatementBlock::addStatement
returnValue addStatement(const ExportStatement &_statement)
Definition: export_statement_block.cpp:91
CG_USE_OPENMP
Definition: acado_types.hpp:339
ExportData::getFullName
std::string getFullName() const
Definition: export_data.cpp:121
REAL
Definition: acado_types.hpp:796
ExportStatementBlock::addLinebreak
returnValue addLinebreak(uint num=1)
Definition: export_statement_block.cpp:209
ExportFunction::setReturnValue
ExportFunction & setReturnValue(const ExportVariable &_functionReturnValue, bool _returnAsPointer=false)
Definition: export_function.cpp:141
ExportGaussElim::getDataDeclarations
virtual returnValue getDataDeclarations(ExportStatementBlock &declarations, ExportStruct dataStruct=ACADO_ANY) const
Definition: gaussian_elimination_export.cpp:53
ExportGaussElim::setupSolveUpperTriangular
virtual returnValue setupSolveUpperTriangular(ExportFunction &_solveTriangular)
Definition: gaussian_elimination_export.cpp:123
ACADO_LOCAL
Definition: acado_types.hpp:808
ExportStatementBlock::addDeclaration
returnValue addDeclaration(const ExportVariable &_data, ExportStruct _dataStruct=ACADO_ANY)
Definition: export_statement_block.cpp:158
ExportGaussElim::setupSolve
virtual returnValue setupSolve(ExportFunction &_solve, ExportFunction &_solveTriangular, ExportVariable &_swap, ExportVariable &_determinant, const std::string &absF)
Definition: gaussian_elimination_export.cpp:141
ExportGaussElim::ExportGaussElim
ExportGaussElim(UserInteraction *_userInteraction=0, const std::string &_commonHeaderName="")
Definition: gaussian_elimination_export.cpp:44
ExportLinearSolver::getNameSolveFunction
const std::string getNameSolveFunction()
Definition: linear_solver_export.cpp:175
BEGIN_NAMESPACE_ACADO
#define BEGIN_NAMESPACE_ACADO
Definition: acado_namespace_macros.hpp:57
ExportStatementBlock::addFunction
returnValue addFunction(const ExportFunction &_function)
Definition: export_statement_block.cpp:108
ExportGaussElim::setupSolveReuse
virtual returnValue setupSolveReuse(ExportFunction &_solveReuse, ExportFunction &_solveTriangular, ExportVariable &_bPerm)
Definition: gaussian_elimination_export.cpp:391
ExportLinearSolver
Allows to export automatically generated algorithms for solving linear systems of specific dimensions...
Definition: linear_solver_export.hpp:62
ExportGaussElim::getFunctionDeclarations
virtual returnValue getFunctionDeclarations(ExportStatementBlock &declarations) const
Definition: gaussian_elimination_export.cpp:69
ExportData::getType
ExportType getType() const
Definition: export_data.cpp:91
ExportStatementBlock
Allows to export code for a block of statements.
Definition: export_statement_block.hpp:63
ExportGaussElim::getGlobalExportVariable
virtual ExportVariable getGlobalExportVariable(const uint factor) const
Definition: gaussian_elimination_export.cpp:577
ExportVariable::getCol
ExportVariable getCol(const ExportIndex &idx) const
Definition: export_variable.cpp:402
ExportFunction::addIndex
ExportFunction & addIndex(const ExportIndex &_index)
Definition: export_function.cpp:347
ACADOERROR
#define ACADOERROR(retval)
Definition: acado_message_handling.hpp:106
ExportVariable
Defines a matrix-valued variable to be used for exporting code.
Definition: export_variable.hpp:58
ExportStatementBlock::addFunctionCall
returnValue addFunctionCall(const std::string &_fName, const ExportArgument &_argument1=emptyConstExportArgument, const ExportArgument &_argument2=emptyConstExportArgument, const ExportArgument &_argument3=emptyConstExportArgument, const ExportArgument &_argument4=emptyConstExportArgument, const ExportArgument &_argument5=emptyConstExportArgument, const ExportArgument &_argument6=emptyConstExportArgument, const ExportArgument &_argument7=emptyConstExportArgument, const ExportArgument &_argument8=emptyConstExportArgument, const ExportArgument &_argument9=emptyConstExportArgument)
Definition: export_statement_block.cpp:115
ExportData::getName
std::string getName() const
Definition: export_data.cpp:86
ExportGaussElim::setupFactorization
virtual returnValue setupFactorization(ExportFunction &_solve, ExportVariable &_swap, ExportVariable &_determinant, const std::string &absF)
Definition: gaussian_elimination_export.cpp:273

acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:37