symmetric_conjugate_gradient_method.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 
35 
36 
38 
39 
40 
41 //
42 // PUBLIC MEMBER FUNCTIONS:
43 //
44 
45 
48 
49  index = 0;
50  nIndex = 0;
51  diag = 0;
52 }
53 
54 
57 
58  int run1, run2;
59 
60  if( arg.nIndex != 0 ){
61  nIndex = new int[dim];
62  for( run1 = 0; run1 < dim; run1++ ){
63  nIndex[run1] = arg.nIndex[run1];
64  }
65  }
66  else nIndex = 0;
67 
68  if( arg.index != 0 ){
69  index = new int*[dim];
70  for( run1 = 0; run1 < dim; run1++ ){
71  index[run1] = new int[nIndex[run1]];
72  for( run2 = 0; run2 < nIndex[run1]; run2++ )
73  index[run1][run2] = arg.index[run1][run2];
74  }
75  }
76  else index = 0;
77 
78  if( arg.diag != 0 ){
79  diag = new int[dim];
80  for( run1 = 0; run1 < dim; run1++ ){
81  diag[run1] = arg.diag[run1];
82  }
83  }
84  else diag = 0;
85 }
86 
87 
89 
90  int run1;
91 
92  if( nIndex != 0 ) delete[] nIndex;
93 
94  if( diag != 0 ) delete[] diag;
95 
96  if( index != 0 ){
97  for( run1 = 0; run1 < dim; run1++ )
98  delete[] index[run1];
99  delete[] index;
100  }
101 }
102 
103 
104 
106 
107  return new SymmetricConjugateGradientMethod(*this);
108 }
109 
110 
111 returnValue SymmetricConjugateGradientMethod::setIndices( const int *rowIdx_, const int *colIdx_ ){
112 
114 }
115 
116 
118 
119 
120  // CONSISTENCY CHECK:
121  // ------------------
122 
123  if( dim <= 0 ) return ACADOERROR(RET_MEMBER_NOT_INITIALISED);
124  if( nDense <= 0 ) return ACADOERROR(RET_MEMBER_NOT_INITIALISED);
125 
126 
127  // ALLOCATE MEMORY:
128  // ------------------
129 
130  if( nIndex != 0 ) delete[] nIndex;
131  nIndex = new int[dim];
132 
133  if( diag != 0 ) delete[] diag;
134  diag = new int[dim];
135 
136 
137  // LOCAL AUXILIARY VARIABLE:
138  // -------------------------
139  int run1, run2;
140  int counter0 ;
141  int bound ;
142  int bound2 ;
143  int counter = 0;
144 
145  BooleanType diagFound;
146 
147 
148  // DETERMINE THE INDEX LENGHTS:
149  // ----------------------------
150 
151 
152  for( run1 = 0; run1 < dim; run1++ ){
153 
154  counter0 = counter;
155  bound = dim*(run1+1);
156  bound2 = run1*(dim+1);
157 
158  diagFound = BT_FALSE;
159 
160  while( counter < nDense ){
161 
162  if( diagFound == BT_FALSE && indices[counter] >= bound2 ){
163  diag[run1] = counter;
164  diagFound = BT_TRUE;
165  }
166 
167  if( indices[counter] >= bound ) break;
168  counter++;
169  }
170 
171  nIndex[run1] = counter - counter0;
172  }
173 
174  if( index != 0 ){
175  for( run1 = 0; run1 < dim; run1++ )
176  delete[] index[run1];
177  delete[] index;
178  }
179  index = new int*[dim];
180 
181  counter = 0;
182 
183  for( run1 = 0; run1 < dim; run1++ ){
184  index[run1] = new int[nIndex[run1]];
185  for( run2 = 0; run2 < nIndex[run1]; run2++ ){
186  index [run1][run2] = indices[counter];
187  counter++;
188  }
189  }
190 
191  return SUCCESSFUL_RETURN;
192 }
193 
194 
195 
196 
197 
198 
199 
200 
201 //
202 // PROTECTED MEMBER FUNCTIONS:
203 //
204 
205 
206 
207 void SymmetricConjugateGradientMethod::multiply( double *xx , double *result ){
208 
209 
210  int i,j,counter, offset;
211 
212  counter = 0;
213 
214  for( i = 0; i < dim; i++ ){
215  result[i] = 0.0;
216  offset = dim*i;
217  for( j = 0; j < nIndex[i]; j++ ){
218  result[i] += A[counter]*xx[index[i][j]-offset];
219  counter++;
220  }
221  }
222 }
223 
224 
225 
226 
228 
229  int i,j,counter, offset;
230 
231  for( i = 0; i < dim; i++ ){
232  condScale[i] = sqrt(A_[diag[i]]);
233  }
234 
235  counter = 0;
236 
237  for( i = 0; i < dim; i++ ){
238  offset = dim*i;
239  for( j = 0; j < nIndex[i]; j++ ){
240  A[counter] = A_[counter]/(condScale[i]*condScale[index[i][j]-offset]);
241  counter++;
242  }
243  }
244 
245  return SUCCESSFUL_RETURN;
246 }
247 
248 
249 
251 
252  int run1;
253 
254  for( run1 = 0; run1 < dim; run1++ )
255  r[run1] = b[run1]/condScale[run1];
256 
257  return SUCCESSFUL_RETURN;
258 }
259 
260 
261 
262 
264 
265  int run1;
266 
267  for( run1 = 0; run1 < dim; run1++ )
268  x_[run1] = x[run1]/condScale[run1];
269 
270  return SUCCESSFUL_RETURN;
271 }
272 
273 
274 
275 
277 
278 
279 /*
280  * end of file
281  */
int counter
Definition: powerkite_c.cpp:40
IntermediateState sqrt(const Expression &arg)
virtual returnValue applyInversePreconditioner(double *x_)
Allows to pass back messages to the calling function.
virtual returnValue setIndices(const int *rowIdx_, const int *colIdx_)
#define CLOSE_NAMESPACE_ACADO
Implements a conjugate gradient method as sparse linear algebra solver for symmetric linear systems...
virtual void multiply(double *xx, double *result)
Implements a conjugate gradient method as sparse linear algebra solver.
#define BT_TRUE
Definition: acado_types.hpp:47
Generic interface for sparse solvers to be coupled with ACADO Toolkit.
#define BEGIN_NAMESPACE_ACADO
#define BT_FALSE
Definition: acado_types.hpp:49
#define ACADOERROR(retval)


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