qpOASES-3.0beta/src/Matrices.cpp
Go to the documentation of this file.
1 /*
2  * This file is part of qpOASES.
3  *
4  * qpOASES -- An Implementation of the Online Active Set Strategy.
5  * Copyright (C) 2007-2011 by Hans Joachim Ferreau, Andreas Potschka,
6  * Christian Kirches et al. All rights reserved.
7  *
8  * qpOASES is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * qpOASES is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with qpOASES; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  *
22  */
23 
24 
35 #include <qpOASES.hpp>
36 
37 #include <qpOASES/PrivateUtils.hpp>
38 
40 
42 const char * const TRANS = "TRANS";
43 
45 const char * const NOTRANS = "NOTRANS";
46 
47 /*****************************************************************************
48  * P U B L I C *
49  *****************************************************************************/
50 
52 {
53  if ( needToFreeMemory( ) == BT_TRUE )
54  free( );
55 }
56 
58 {
59  if (val != 0)
60  delete[] val;
61  val = 0;
62 }
63 
65 {
66  DenseMatrix *dupl = 0;
67 
68  if ( needToFreeMemory( ) == BT_TRUE )
69  {
70  real_t* val_new = new real_t[nRows*nCols];
71  memcpy( val_new,val, nRows*nCols*sizeof(real_t) );
72  dupl = new DenseMatrix(nRows, nCols, nCols, val_new);
73  dupl->doFreeMemory( );
74  }
75  else
76  {
77  dupl = new DenseMatrix(nRows, nCols, nCols, val);
78  }
79 
80  return dupl;
81 }
82 
84  ) const
85 {
86  return val[i*(leaDim+1)];
87 }
88 
90 {
91  int i, j;
92 
93  if (nRows != nCols)
94  return BT_FALSE;
95 
96  for ( i=0; i<nRows; ++i )
97  for ( j=0; j<i; ++j )
98  if ( ( fabs( val[i*leaDim+j] ) > EPS ) || ( fabs( val[j*leaDim+i] ) > EPS ) )
99  return BT_FALSE;
100 
101  return BT_TRUE;
102 }
103 
104 returnValue DenseMatrix::getRow(int rNum, const Indexlist* const icols, real_t alpha, real_t *row) const
105 {
106  int i;
107  if (icols != 0)
108  {
109  if (isExactlyOne(alpha))
110  for (i = 0; i < icols->length; i++)
111  row[i] = val[rNum*leaDim+icols->number[i]];
112  else if (isExactlyMinusOne(alpha))
113  for (i = 0; i < icols->length; i++)
114  row[i] = -val[rNum*leaDim+icols->number[i]];
115  else
116  for (i = 0; i < icols->length; i++)
117  row[i] = alpha*val[rNum*leaDim+icols->number[i]];
118  }
119  else
120  {
121  if (isExactlyOne(alpha))
122  for (i = 0; i < nCols; i++)
123  row[i] = val[rNum*leaDim+i];
124  else if (isExactlyMinusOne(alpha))
125  for (i = 0; i < nCols; i++)
126  row[i] = -val[rNum*leaDim+i];
127  else
128  for (i = 0; i < nCols; i++)
129  row[i] = alpha*val[rNum*leaDim+i];
130  }
131  return SUCCESSFUL_RETURN;
132 }
133 
134 returnValue DenseMatrix::getCol(int cNum, const Indexlist* const irows, real_t alpha, real_t *col) const
135 {
136  int i;
137 
138  if (isExactlyOne(alpha))
139  for (i = 0; i < irows->length; i++)
140  col[i] = val[irows->number[i]*leaDim+cNum];
141  else if (isExactlyMinusOne(alpha))
142  for (i = 0; i < irows->length; i++)
143  col[i] = -val[irows->number[i]*leaDim+cNum];
144  else
145  for (i = 0; i < irows->length; i++)
146  col[i] = alpha*val[irows->number[i]*leaDim+cNum];
147 
148  return SUCCESSFUL_RETURN;
149 }
150 
151 returnValue DenseMatrix::times(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
152 {
153  unsigned long _xN=xN, _nRows=nRows, _nCols=nCols, /*_leaDim=getMax(1,leaDim),*/ _xLD=getMax(1,xLD), _yLD=getMax(1,yLD);
154  /* Call BLAS. Mind row major format! */
155  GEMM(TRANS, NOTRANS, &_nRows, &_xN, &_nCols, &alpha, val, &_nCols, x, &_xLD, &beta, y, &_yLD);
156  return SUCCESSFUL_RETURN;
157 }
158 
159 returnValue DenseMatrix::transTimes(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
160 {
161  /* Call BLAS. Mind row major format! */
162  unsigned long _xN=xN, _nRows=nRows, _nCols=nCols, /*_leaDim=getMax(1,leaDim),*/ _xLD=getMax(1,xLD), _yLD=getMax(1,yLD);
163  GEMM(NOTRANS, NOTRANS, &_nCols, &_xN, &_nRows, &alpha, val, &_nCols, x, &_xLD, &beta, y, &_yLD);
164  return SUCCESSFUL_RETURN;
165 }
166 
167 returnValue DenseMatrix::times(const Indexlist* const irows, const Indexlist* const icols,
168  int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD,
169  BooleanType yCompr) const
170 {
171  int i, j, k, row, col, iy, irA;
172 
173  if (yCompr == BT_TRUE)
174  {
175  if (isExactlyZero(beta))
176  for (k = 0; k < xN; k++)
177  for (j = 0; j < irows->length; j++)
178  y[j+k*yLD] = 0.0;
179  else if (isExactlyMinusOne(beta))
180  for (k = 0; k < xN; k++)
181  for (j = 0; j < irows->length; j++)
182  y[j+k*yLD] = -y[j+k*yLD];
183  else if (!isExactlyOne(beta))
184  for (k = 0; k < xN; k++)
185  for (j = 0; j < irows->length; j++)
186  y[j+k*yLD] *= beta;
187 
188  if (icols == 0)
189  if (isExactlyOne(alpha))
190  for (k = 0; k < xN; k++)
191  for (j = 0; j < irows->length; j++)
192  {
193  row = irows->iSort[j];
194  iy = row + k * yLD;
195  irA = irows->number[row] * leaDim;
196  for (i = 0; i < nCols; i++)
197  y[iy] += val[irA+i] * x[k*xLD+i];
198  }
199  else if (isExactlyMinusOne(alpha))
200  for (k = 0; k < xN; k++)
201  for (j = 0; j < irows->length; j++)
202  {
203  row = irows->iSort[j];
204  iy = row + k * yLD;
205  irA = irows->number[row] * leaDim;
206  for (i = 0; i < nCols; i++)
207  y[iy] -= val[irA+i] * x[k*xLD+i];
208  }
209  else
210  for (k = 0; k < xN; k++)
211  for (j = 0; j < irows->length; j++)
212  {
213  row = irows->iSort[j];
214  iy = row + k * yLD;
215  irA = irows->number[row] * leaDim;
216  for (i = 0; i < nCols; i++)
217  y[iy] += alpha * val[irA+i] * x[k*xLD+i];
218  }
219  else /* icols != 0 */
220  if (isExactlyOne(alpha))
221  for (k = 0; k < xN; k++)
222  for (j = 0; j < irows->length; j++)
223  {
224  row = irows->iSort[j];
225  iy = row + k * yLD;
226  irA = irows->number[row] * leaDim;
227  for (i = 0; i < icols->length; i++)
228  {
229  col = icols->iSort[i];
230  y[iy] += val[irA+icols->number[col]] * x[k*xLD+col];
231  }
232  }
233  else if (isExactlyMinusOne(alpha))
234  for (k = 0; k < xN; k++)
235  for (j = 0; j < irows->length; j++)
236  {
237  row = irows->iSort[j];
238  iy = row + k * yLD;
239  irA = irows->number[row] * leaDim;
240  for (i = 0; i < icols->length; i++)
241  {
242  col = icols->iSort[i];
243  y[iy] -= val[irA+icols->number[col]] * x[k*xLD+col];
244  }
245  }
246  else
247  for (k = 0; k < xN; k++)
248  for (j = 0; j < irows->length; j++)
249  {
250  row = irows->iSort[j];
251  iy = row + k * yLD;
252  irA = irows->number[row] * leaDim;
253  for (i = 0; i < icols->length; i++)
254  {
255  col = icols->iSort[i];
256  y[iy] += alpha * val[irA+icols->number[col]] * x[k*xLD+col];
257  }
258  }
259  }
260  else /* y not compressed */
261  {
262  if (isExactlyZero(beta))
263  for (k = 0; k < xN; k++)
264  for (j = 0; j < irows->length; j++)
265  y[irows->number[j]+k*yLD] = 0.0;
266  else if (isExactlyMinusOne(beta))
267  for (k = 0; k < xN; k++)
268  for (j = 0; j < irows->length; j++)
269  y[irows->number[j]+k*yLD] = -y[j+k*yLD];
270  else if (!isExactlyOne(beta))
271  for (k = 0; k < xN; k++)
272  for (j = 0; j < irows->length; j++)
273  y[irows->number[j]+k*yLD] *= beta;
274 
275  if (icols == 0)
276  if (isExactlyOne(alpha))
277  for (k = 0; k < xN; k++)
278  for (j = 0; j < irows->length; j++)
279  {
280  row = irows->number[irows->iSort[j]];
281  iy = row + k * yLD;
282  irA = row * leaDim;
283  for (i = 0; i < nCols; i++)
284  y[iy] += val[irA+i] * x[k*xLD+i];
285  }
286  else if (isExactlyMinusOne(alpha))
287  for (k = 0; k < xN; k++)
288  for (j = 0; j < irows->length; j++)
289  {
290  row = irows->number[irows->iSort[j]];
291  iy = row + k * yLD;
292  irA = row * leaDim;
293  for (i = 0; i < nCols; i++)
294  y[iy] -= val[irA+i] * x[k*xLD+i];
295  }
296  else
297  for (k = 0; k < xN; k++)
298  for (j = 0; j < irows->length; j++)
299  {
300  row = irows->number[irows->iSort[j]];
301  iy = row + k * yLD;
302  irA = row * leaDim;
303  for (i = 0; i < nCols; i++)
304  y[iy] += alpha * val[irA+i] * x[k*xLD+i];
305  }
306  else /* icols != 0 */
307  if (isExactlyOne(alpha))
308  for (k = 0; k < xN; k++)
309  for (j = 0; j < irows->length; j++)
310  {
311  row = irows->number[irows->iSort[j]];
312  iy = row + k * yLD;
313  irA = row * leaDim;
314  for (i = 0; i < icols->length; i++)
315  {
316  col = icols->iSort[i];
317  y[iy] += val[irA+icols->number[col]] * x[k*xLD+col];
318  }
319  }
320  else if (isExactlyMinusOne(alpha))
321  for (k = 0; k < xN; k++)
322  for (j = 0; j < irows->length; j++)
323  {
324  row = irows->number[irows->iSort[j]];
325  iy = row + k * yLD;
326  irA = row * leaDim;
327  for (i = 0; i < icols->length; i++)
328  {
329  col = icols->iSort[i];
330  y[iy] -= val[irA+icols->number[col]] * x[k*xLD+col];
331  }
332  }
333  else
334  for (k = 0; k < xN; k++)
335  for (j = 0; j < irows->length; j++)
336  {
337  row = irows->number[irows->iSort[j]];
338  iy = row + k * yLD;
339  irA = row * leaDim;
340  for (i = 0; i < icols->length; i++)
341  {
342  col = icols->iSort[i];
343  y[iy] += alpha * val[irA+icols->number[col]] * x[k*xLD+col];
344  }
345  }
346  }
347 
348  return SUCCESSFUL_RETURN;
349 }
350 
351 returnValue DenseMatrix::transTimes(const Indexlist* const irows, const Indexlist* const icols,
352  int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
353 {
354  int i, j, k, row, col;
355 
356  if (isExactlyZero(beta))
357  for (k = 0; k < xN; k++)
358  for (j = 0; j < icols->length; j++)
359  y[j+k*yLD] = 0.0;
360  else if (isExactlyMinusOne(beta))
361  for (k = 0; k < xN; k++)
362  for (j = 0; j < icols->length; j++)
363  y[j+k*yLD] = -y[j+k*yLD];
364  else if (!isExactlyOne(beta))
365  for (k = 0; k < xN; k++)
366  for (j = 0; j < icols->length; j++)
367  y[j+k*yLD] *= beta;
368 
369  if (isExactlyOne(alpha))
370  for (k = 0; k < xN; k++)
371  for (j = 0; j < irows->length; j++)
372  {
373  row = irows->iSort[j];
374  for (i = 0; i < icols->length; i++)
375  {
376  col = icols->iSort[i];
377  y[col+k*yLD] += val[irows->number[row]*leaDim+icols->number[col]] * x[row+k*xLD];
378  }
379  }
380  else if (isExactlyMinusOne(alpha))
381  for (k = 0; k < xN; k++)
382  for (j = 0; j < irows->length; j++)
383  {
384  row = irows->iSort[j];
385  for (i = 0; i < icols->length; i++)
386  {
387  col = icols->iSort[i];
388  y[col+k*yLD] -= val[irows->number[row]*leaDim+icols->number[col]] * x[row+k*xLD];
389  }
390  }
391  else
392  for (k = 0; k < xN; k++)
393  for (j = 0; j < irows->length; j++)
394  {
395  row = irows->iSort[j];
396  for (i = 0; i < icols->length; i++)
397  {
398  col = icols->iSort[i];
399  y[col+k*yLD] += alpha * val[irows->number[row]*leaDim+icols->number[col]] * x[row+k*xLD];
400  }
401  }
402 
403  return SUCCESSFUL_RETURN;
404 }
405 
407 {
408  int i;
409  for (i = 0; i < nRows && i < nCols; i++)
410  val[i*(leaDim+1)] += alpha;
411 
412  return SUCCESSFUL_RETURN;
413 }
414 
415 
417 {
418  myPrintf( "aaa\n" );
419  return qpOASES::print( val,nRows,nCols );
420 }
421 
422 
423 
425 {
426  /* "same" as duplicate() in DenseMatrix */
427  SymDenseMat *dupl = 0;
428 
429  if ( needToFreeMemory( ) == BT_TRUE )
430  {
431  real_t* val_new = new real_t[nRows*nCols];
432  memcpy( val_new,val, nRows*nCols*sizeof(real_t) );
433  dupl = new SymDenseMat(nRows, nCols, nCols, val_new);
434  dupl->doFreeMemory( );
435  }
436  else
437  {
438  dupl = new SymDenseMat(nRows, nCols, nCols, val);
439  }
440 
441  return dupl;
442 }
443 
444 
446  int xN, const real_t *x, int xLD, real_t *y, int yLD) const
447 {
448  int ii, jj, kk, col;
449  int i,j,k,irA;
450 
451  for (ii = 0; ii < xN; ii++)
452  for (jj = 0; jj < xN; jj++)
453  y[ii*yLD+jj] = 0.0;
454 
455  real_t *Ax = new real_t[icols->length * xN];
456 
457  for (i=0;i<icols->length * xN;++i)
458  Ax[i]=0.0;
459 
460  /* exploit symmetry of A ! */
461  for (j = 0; j < icols->length; j++) {
462  irA = icols->number[j] * leaDim;
463  for (i = 0; i < icols->length; i++)
464  {
465  real_t h = val[irA+icols->number[i]];
466  for (k = 0; k < xN; k++)
467  Ax[j + k * icols->length] += h * x[k*xLD+icols->number[i]];
468  }
469  }
470 
471  for (ii = 0; ii < icols->length; ++ii) {
472  col = icols->number[ii];
473  for (jj = 0; jj < xN; ++jj) {
474  for (kk = 0; kk < xN; ++kk) {
475  y[kk + jj*yLD] += x[col + jj*xLD] * Ax[ii + kk*icols->length];
476  }
477  }
478  }
479  delete[] Ax;
480 
481  return SUCCESSFUL_RETURN;
482 }
483 
484 
485 SparseMatrix::SparseMatrix() : nRows(0), nCols(0), ir(0), jc(0), jd(0), val(0) {}
486 
487 SparseMatrix::SparseMatrix(long nr, long nc, long *r, long *c, real_t *v, long *d)
488  : nRows(nr), nCols(nc), ir(r), jc(c), jd(d), val(v) { doNotFreeMemory(); }
489 
490 SparseMatrix::SparseMatrix(int nr, int nc, int ld, const real_t * const v) : nRows(nr), nCols(nc), jd(0)
491 {
492  int i, j, nnz;
493 
494  jc = new long[nc+1];
495  ir = new long[nr*nc];
496  val = new real_t[nr*nc];
497 
498  nnz = 0;
499  for (j = 0; j < nCols; j++)
500  {
501  jc[j] = nnz;
502  for (i = 0; i < nRows; i++)
503  if (!isExactlyZero(v[i*ld+j]))
504  {
505  ir[nnz] = i;
506  val[nnz++] = v[i*ld+j];
507  }
508  }
509  jc[nCols] = nnz;
510 
511  doFreeMemory( );
512 }
513 
515 {
516  if ( needToFreeMemory() == BT_TRUE )
517  free( );
518 }
519 
521 {
522  if (ir != 0) delete[] ir;
523  ir = 0;
524  if (jc != 0) delete[] jc;
525  jc = 0;
526  if (jd != 0) delete[] jd;
527  jd = 0;
528  if (val != 0) delete[] val;
529  val = 0;
530  doNotFreeMemory( );
531 }
532 
534 {
535  long i, length = jc[nCols];
536  SparseMatrix *dupl = new SparseMatrix;
537 
538  dupl->nRows = nRows;
539  dupl->nCols = nCols;
540  dupl->ir = new long[length];
541  dupl->jc = new long[nCols+1];
542  dupl->jd = new long[nCols];
543  dupl->val = new real_t[length];
544 
545  for (i = 0; i < length; i++) dupl->ir[i] = ir[i];
546  for (i = 0; i <= nCols; i++) dupl->jc[i] = jc[i];
547  for (i = 0; i < nCols; i++) dupl->jd[i] = jd[i];
548  for (i = 0; i < length; i++) dupl->val[i] = val[i];
549 
550  dupl->doFreeMemory( );
551 
552  return dupl;
553 }
554 
556 {
557  int entry = jd[i];
558  return (entry < jc[i+1] && ir[entry] == i) ? val[entry] : 0.0;
559 }
560 
562 {
563  return BT_FALSE;
564 }
565 
566 returnValue SparseMatrix::getRow(int rNum, const Indexlist* const icols, real_t alpha, real_t *row) const
567 {
568  long i, j, k;
569 
570  if (icols != 0)
571  {
572  if (isExactlyOne(alpha))
573  for (k = 0; k < icols->length; k++)
574  {
575  j = icols->number[icols->iSort[k]];
576  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
577  row[icols->iSort[k]] = (i < jc[j+1] && ir[i] == rNum) ? val[i] : 0.0;
578  }
579  else if (isExactlyMinusOne(alpha))
580  for (k = 0; k < icols->length; k++)
581  {
582  j = icols->number[icols->iSort[k]];
583  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
584  row[icols->iSort[k]] = (i < jc[j+1] && ir[i] == rNum) ? -val[i] : 0.0;
585  }
586  else
587  for (k = 0; k < icols->length; k++)
588  {
589  j = icols->number[icols->iSort[k]];
590  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
591  row[icols->iSort[k]] = (i < jc[j+1] && ir[i] == rNum) ? alpha*val[i] : 0.0;
592  }
593  }
594  else
595  {
596  if (isExactlyOne(alpha))
597  for (j = 0; j < nCols; j++)
598  {
599  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
600  row[j] = (i < jc[j+1] && ir[i] == rNum) ? val[i] : 0.0;
601  }
602  else if (isExactlyMinusOne(alpha))
603  for (j = 0; j < icols->length; j++)
604  {
605  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
606  row[j] = (i < jc[j+1] && ir[i] == rNum) ? -val[i] : 0.0;
607  }
608  else
609  for (j = 0; j < icols->length; j++)
610  {
611  for (i = jc[j]; i < jc[j+1] && ir[i] < rNum; i++);
612  row[j] = (i < jc[j+1] && ir[i] == rNum) ? alpha*val[i] : 0.0;
613  }
614  }
615  return SUCCESSFUL_RETURN;
616 }
617 
618 returnValue SparseMatrix::getCol(int cNum, const Indexlist* const irows, real_t alpha, real_t *col) const
619 {
620  long i, j;
621 
622  i = jc[cNum];
623  j = 0;
624  if (isExactlyOne(alpha))
625  while (i < jc[cNum+1] && j < irows->length)
626  if (ir[i] == irows->number[irows->iSort[j]])
627  col[irows->iSort[j++]] = val[i++];
628  else if (ir[i] > irows->number[irows->iSort[j]])
629  col[irows->iSort[j++]] = 0.0;
630  else
631  i++;
632  else if (isExactlyMinusOne(alpha))
633  while (i < jc[cNum+1] && j < irows->length)
634  if (ir[i] == irows->number[irows->iSort[j]])
635  col[irows->iSort[j++]] = -val[i++];
636  else if (ir[i] > irows->number[irows->iSort[j]])
637  col[irows->iSort[j++]] = 0.0;
638  else
639  i++;
640  else
641  while (i < jc[cNum+1] && j < irows->length)
642  if (ir[i] == irows->number[irows->iSort[j]])
643  col[irows->iSort[j++]] = alpha * val[i++];
644  else if (ir[i] > irows->number[irows->iSort[j]])
645  col[irows->iSort[j++]] = 0.0;
646  else
647  i++;
648 
649  /* fill in remaining zeros */
650  while (j < irows->length)
651  col[irows->iSort[j++]] = 0.0;
652 
653  return SUCCESSFUL_RETURN;
654 }
655 
656 returnValue SparseMatrix::times(int xN, real_t alpha, const real_t *x, int xLD,
657  real_t beta, real_t *y, int yLD) const
658 {
659  long i, j, k;
660 
661  if (isExactlyZero(beta))
662  for (k = 0; k < xN; k++)
663  for (j = 0; j < nRows; j++)
664  y[j+k*yLD] = 0.0;
665  else if (isExactlyMinusOne(beta))
666  for (k = 0; k < xN; k++)
667  for (j = 0; j < nRows; j++)
668  y[j+k*yLD] = -y[j+k*yLD];
669  else if (!isExactlyOne(beta))
670  for (k = 0; k < xN; k++)
671  for (j = 0; j < nRows; j++)
672  y[j+k*yLD] *= beta;
673 
674  if (isExactlyOne(alpha))
675  for (k = 0; k < xN; k++)
676  for (j = 0; j < nCols; j++)
677  for (i = jc[j]; i < jc[j+1]; i++)
678  y[ir[i]+k*yLD] += val[i] * x[j+k*xLD];
679  else if (isExactlyMinusOne(alpha))
680  for (k = 0; k < xN; k++)
681  for (j = 0; j < nCols; j++)
682  for (i = jc[j]; i < jc[j+1]; i++)
683  y[ir[i]+k*yLD] -= val[i] * x[j+k*xLD];
684  else
685  for (k = 0; k < xN; k++)
686  for (j = 0; j < nCols; j++)
687  for (i = jc[j]; i < jc[j+1]; i++)
688  y[ir[i]+k*yLD] += alpha * val[i] * x[j+k*xLD];
689 
690  return SUCCESSFUL_RETURN;
691 }
692 
693 returnValue SparseMatrix::transTimes(int xN, real_t alpha, const real_t *x, int xLD,
694  real_t beta, real_t *y, int yLD) const
695 {
696  long i, j, k;
697 
698  if (isExactlyZero(beta))
699  for (k = 0; k < xN; k++)
700  for (j = 0; j < nCols; j++)
701  y[j+k*yLD] = 0.0;
702  else if (isExactlyMinusOne(beta))
703  for (k = 0; k < xN; k++)
704  for (j = 0; j < nCols; j++)
705  y[j+k*yLD] = -y[j+k*yLD];
706  else if (!isExactlyOne(beta))
707  for (k = 0; k < xN; k++)
708  for (j = 0; j < nCols; j++)
709  y[j+k*yLD] *= beta;
710 
711  if (isExactlyOne(alpha))
712  for (k = 0; k < xN; k++)
713  for (j = 0; j < nCols; j++)
714  for (i = jc[j]; i < jc[j+1]; i++)
715  y[j+k*yLD] += val[i] * x[ir[i]+k*xLD];
716  else if (isExactlyMinusOne(alpha))
717  for (k = 0; k < xN; k++)
718  for (j = 0; j < nCols; j++)
719  for (i = jc[j]; i < jc[j+1]; i++)
720  y[j+k*yLD] -= val[i] * x[ir[i]+k*xLD];
721  else
722  for (k = 0; k < xN; k++)
723  for (j = 0; j < nCols; j++)
724  for (i = jc[j]; i < jc[j+1]; i++)
725  y[j+k*yLD] += alpha * val[i] * x[ir[i]+k*xLD];
726 
727  return SUCCESSFUL_RETURN;
728 }
729 
730 returnValue SparseMatrix::times(const Indexlist* const irows, const Indexlist* const icols,
731  int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD,
732  BooleanType yCompr) const
733 {
734  long i, j, k, l, col;
735 
736  if (yCompr == BT_TRUE)
737  {
738  if (isExactlyZero(beta))
739  for (k = 0; k < xN; k++)
740  for (j = 0; j < irows->length; j++)
741  y[j+k*yLD] = 0.0;
742  else if (isExactlyMinusOne(beta))
743  for (k = 0; k < xN; k++)
744  for (j = 0; j < irows->length; j++)
745  y[j+k*yLD] = -y[j+k*yLD];
746  else if (!isExactlyOne(beta))
747  for (k = 0; k < xN; k++)
748  for (j = 0; j < irows->length; j++)
749  y[j+k*yLD] *= beta;
750 
751  if (icols == 0)
752  if (isExactlyOne(alpha))
753  for (l = 0; l < nCols; l++)
754  {
755  i = jc[l];
756  j = 0;
757  while (i < jc[l+1] && j < irows->length)
758  if (ir[i] == irows->number[irows->iSort[j]])
759  {
760  for (k = 0; k < xN; k++)
761  y[k*yLD+irows->iSort[j]] += val[i] * x[k*xLD+l];
762  i++, j++;
763  }
764  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
765  else i++;
766  }
767  else if (isExactlyMinusOne(alpha))
768  for (l = 0; l < nCols; l++)
769  {
770  i = jc[l];
771  j = 0;
772  while (i < jc[l+1] && j < irows->length)
773  if (ir[i] == irows->number[irows->iSort[j]])
774  {
775  for (k = 0; k < xN; k++)
776  y[k*yLD+irows->iSort[j]] -= val[i] * x[k*xLD+l];
777  i++, j++;
778  }
779  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
780  else i++;
781  }
782  else
783  for (l = 0; l < nCols; l++)
784  {
785  i = jc[l];
786  j = 0;
787  while (i < jc[l+1] && j < irows->length)
788  if (ir[i] == irows->number[irows->iSort[j]])
789  {
790  for (k = 0; k < xN; k++)
791  y[k*yLD+irows->iSort[j]] += alpha * val[i] * x[k*xLD+l];
792  i++, j++;
793  }
794  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
795  else i++;
796  }
797  else /* icols != 0 */
798  if (isExactlyOne(alpha))
799  for (l = 0; l < icols->length; l++)
800  {
801  col = icols->iSort[l];
802  i = jc[icols->number[col]];
803  j = 0;
804  while (i < jc[icols->number[col]+1] && j < irows->length)
805  if (ir[i] == irows->number[irows->iSort[j]])
806  {
807  for (k = 0; k < xN; k++)
808  y[k*yLD+irows->iSort[j]] += val[i] * x[k*xLD+col];
809  i++, j++;
810  }
811  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
812  else i++;
813  }
814  else if (isExactlyMinusOne(alpha))
815  for (l = 0; l < icols->length; l++)
816  {
817  col = icols->iSort[l];
818  i = jc[icols->number[col]];
819  j = 0;
820  while (i < jc[icols->number[col]+1] && j < irows->length)
821  if (ir[i] == irows->number[irows->iSort[j]])
822  {
823  for (k = 0; k < xN; k++)
824  y[k*yLD+irows->iSort[j]] -= val[i] * x[k*xLD+col];
825  i++, j++;
826  }
827  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
828  else i++;
829  }
830  else
831  for (l = 0; l < icols->length; l++)
832  {
833  col = icols->iSort[l];
834  i = jc[icols->number[col]];
835  j = 0;
836  while (i < jc[icols->number[col]+1] && j < irows->length)
837  if (ir[i] == irows->number[irows->iSort[j]])
838  {
839  for (k = 0; k < xN; k++)
840  y[k*yLD+irows->iSort[j]] += alpha * val[i] * x[k*xLD+col];
841  i++, j++;
842  }
843  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
844  else i++;
845  }
846  }
847  else /* y not compressed */
848  {
849  if (isExactlyZero(beta))
850  for (k = 0; k < xN; k++)
851  for (j = 0; j < irows->length; j++)
852  y[irows->number[j]+k*yLD] = 0.0;
853  else if (isExactlyMinusOne(beta))
854  for (k = 0; k < xN; k++)
855  for (j = 0; j < irows->length; j++)
856  y[irows->number[j]+k*yLD] = -y[j+k*yLD];
857  else if (!isExactlyOne(beta))
858  for (k = 0; k < xN; k++)
859  for (j = 0; j < irows->length; j++)
860  y[irows->number[j]+k*yLD] *= beta;
861 
862  if (icols == 0)
863  if (isExactlyOne(alpha))
864  for (l = 0; l < nCols; l++)
865  {
866  i = jc[l];
867  j = 0;
868  while (i < jc[l+1] && j < irows->length)
869  if (ir[i] == irows->number[irows->iSort[j]])
870  {
871  for (k = 0; k < xN; k++)
872  y[k*yLD+irows->number[irows->iSort[j]]] += val[i] * x[k*xLD+l];
873  i++, j++;
874  }
875  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
876  else i++;
877  }
878  else if (isExactlyMinusOne(alpha))
879  for (l = 0; l < nCols; l++)
880  {
881  i = jc[l];
882  j = 0;
883  while (i < jc[l+1] && j < irows->length)
884  if (ir[i] == irows->number[irows->iSort[j]])
885  {
886  for (k = 0; k < xN; k++)
887  y[k*yLD+irows->number[irows->iSort[j]]] -= val[i] * x[k*xLD+l];
888  i++, j++;
889  }
890  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
891  else i++;
892  }
893  else
894  for (l = 0; l < nCols; l++)
895  {
896  i = jc[l];
897  j = 0;
898  while (i < jc[l+1] && j < irows->length)
899  if (ir[i] == irows->number[irows->iSort[j]])
900  {
901  for (k = 0; k < xN; k++)
902  y[k*yLD+irows->number[irows->iSort[j]]] += alpha * val[i] * x[k*xLD+l];
903  i++, j++;
904  }
905  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
906  else i++;
907  }
908  else /* icols != 0 */
909  if (isExactlyOne(alpha))
910  for (l = 0; l < icols->length; l++)
911  {
912  col = icols->iSort[l];
913  i = jc[icols->number[col]];
914  j = 0;
915  while (i < jc[icols->number[col]+1] && j < irows->length)
916  if (ir[i] == irows->number[irows->iSort[j]])
917  {
918  for (k = 0; k < xN; k++)
919  y[k*yLD+irows->number[irows->iSort[j]]] += val[i] * x[k*xLD+col];
920  i++, j++;
921  }
922  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
923  else i++;
924  }
925  else if (isExactlyMinusOne(alpha))
926  for (l = 0; l < icols->length; l++)
927  {
928  col = icols->iSort[l];
929  i = jc[icols->number[col]];
930  j = 0;
931  while (i < jc[icols->number[col]+1] && j < irows->length)
932  if (ir[i] == irows->number[irows->iSort[j]])
933  {
934  for (k = 0; k < xN; k++)
935  y[k*yLD+irows->number[irows->iSort[j]]] -= val[i] * x[k*xLD+col];
936  i++, j++;
937  }
938  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
939  else i++;
940  }
941  else
942  for (l = 0; l < icols->length; l++)
943  {
944  col = icols->iSort[l];
945  i = jc[icols->number[col]];
946  j = 0;
947  while (i < jc[icols->number[col]+1] && j < irows->length)
948  if (ir[i] == irows->number[irows->iSort[j]])
949  {
950  for (k = 0; k < xN; k++)
951  y[k*yLD+irows->number[irows->iSort[j]]] += alpha * val[i] * x[k*xLD+col];
952  i++, j++;
953  }
954  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
955  else i++;
956  }
957  }
958  return SUCCESSFUL_RETURN;
959 }
960 
961 returnValue SparseMatrix::transTimes(const Indexlist* const irows, const Indexlist* const icols,
962  int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
963 {
964  long i, j, k, l, col;
965 
966  if (isExactlyZero(beta))
967  for (k = 0; k < xN; k++)
968  for (j = 0; j < icols->length; j++)
969  y[j+k*yLD] = 0.0;
970  else if (isExactlyMinusOne(beta))
971  for (k = 0; k < xN; k++)
972  for (j = 0; j < icols->length; j++)
973  y[j+k*yLD] = -y[j+k*yLD];
974  else if (!isExactlyOne(beta))
975  for (k = 0; k < xN; k++)
976  for (j = 0; j < icols->length; j++)
977  y[j+k*yLD] *= beta;
978 
979  if (isExactlyOne(alpha))
980  for (l = 0; l < icols->length; l++)
981  {
982  col = icols->iSort[l];
983  i = jc[icols->number[col]];
984  j = 0;
985  while (i < jc[icols->number[col]+1] && j < irows->length)
986  if (ir[i] == irows->number[irows->iSort[j]])
987  {
988  for (k = 0; k < xN; k++)
989  y[k*yLD+col] += val[i] * x[k*xLD+irows->iSort[j]];
990  i++, j++;
991  }
992  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
993  else i++;
994  }
995  else if (isExactlyMinusOne(alpha))
996  for (l = 0; l < icols->length; l++)
997  {
998  col = icols->iSort[l];
999  i = jc[icols->number[col]];
1000  j = 0;
1001  while (i < jc[icols->number[col]+1] && j < irows->length)
1002  if (ir[i] == irows->number[irows->iSort[j]])
1003  {
1004  for (k = 0; k < xN; k++)
1005  y[k*yLD+col] -= val[i] * x[k*xLD+irows->iSort[j]];
1006  i++, j++;
1007  }
1008  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
1009  else i++;
1010  }
1011  else
1012  for (l = 0; l < icols->length; l++)
1013  {
1014  col = icols->iSort[l];
1015  i = jc[icols->number[col]];
1016  j = 0;
1017  while (i < jc[icols->number[col]+1] && j < irows->length)
1018  if (ir[i] == irows->number[irows->iSort[j]])
1019  {
1020  for (k = 0; k < xN; k++)
1021  y[k*yLD+col] += alpha * val[i] * x[k*xLD+irows->iSort[j]];
1022  i++, j++;
1023  }
1024  else if (ir[i] > irows->number[irows->iSort[j]]) j++;
1025  else i++;
1026  }
1027 
1028  return SUCCESSFUL_RETURN;
1029 }
1030 
1032 {
1033  long i;
1034 
1035  if (!isExactlyZero(alpha))
1036  {
1037  for (i = 0; i < nRows && i < nCols; i++)
1038  if (ir[jd[i]] == i) val[jd[i]] += alpha;
1039  else return RET_NO_DIAGONAL_AVAILABLE;
1040  }
1041 
1042  return SUCCESSFUL_RETURN;
1043 }
1044 
1046 {
1047  long i, j;
1048 
1049  if (jd == 0) {
1050  jd = new long[nCols];
1051 
1052  for (j = 0; j < nCols; j++)
1053  {
1054  for (i = jc[j]; i < jc[j+1] && ir[i] < j; i++);
1055  jd[j] = i;
1056  }
1057  }
1058 
1059  return jd;
1060 }
1061 
1063 {
1064  long i, j;
1065  real_t *v = new real_t[nRows*nCols];
1066 
1067  for (i = 0; i < nCols*nRows; i++)
1068  v[i] = 0.0;
1069 
1070  for (j = 0; j < nCols; j++)
1071  for (i = jc[j]; i < jc[j+1]; i++)
1072  v[ir[i] * nCols + j] = val[i];
1073 
1074  return v;
1075 }
1076 
1077 
1079 {
1081 }
1082 
1083 
1084 
1086 {
1087  /* "same" as duplicate() in SparseMatrix */
1088  long i, length = jc[nCols];
1089  SymSparseMat *dupl = new SymSparseMat;
1090 
1091  dupl->nRows = nRows;
1092  dupl->nCols = nCols;
1093  dupl->ir = new long[length];
1094  dupl->jc = new long[nCols+1];
1095  dupl->jd = new long[nCols];
1096  dupl->val = new real_t[length];
1097 
1098  for (i = 0; i < length; i++) dupl->ir[i] = ir[i];
1099  for (i = 0; i <= nCols; i++) dupl->jc[i] = jc[i];
1100  for (i = 0; i < nCols; i++) dupl->jd[i] = jd[i];
1101  for (i = 0; i < length; i++) dupl->val[i] = val[i];
1102 
1103  dupl->doFreeMemory( );
1104 
1105  return dupl;
1106 }
1107 
1109  int xN, const real_t *x, int xLD, real_t *y, int yLD) const
1110 {
1111  int i, j, k, l, idx, row, col;
1112 
1113  /* clear output */
1114  for (i = 0; i < xN*xN; i++)
1115  y[i] = 0.0;
1116 
1117  /* compute lower triangle */
1118  for (l = 0; l < icols->length; l++)
1119  {
1120  col = icols->number[icols->iSort[l]];
1121  idx = jd[col];
1122  k = 0;
1123  while (idx < jc[col+1] && k < icols->length)
1124  {
1125  row = icols->number[icols->iSort[k]];
1126  if (ir[idx] == row)
1127  {
1128  /* TODO: It is possible to formulate this as DSYR and DSYR2
1129  * operations. */
1130  if (row == col) /* diagonal element */
1131  for (i = 0; i < xN; i++)
1132  for (j = i; j < xN; j++)
1133  y[i*yLD+j] += val[idx] * x[i*xLD+col] * x[j*xLD+col];
1134  else /* subdiagonal elements */
1135  for (i = 0; i < xN; i++)
1136  for (j = i; j < xN; j++)
1137  y[i*yLD+j] += val[idx] * (x[i*xLD+col] * x[j*xLD+row] + x[i*xLD+row] * x[j*xLD+col]);
1138  idx++, k++;
1139  }
1140  else if (ir[idx] > row) k++;
1141  else idx++;
1142  }
1143  }
1144 
1145  /* fill upper triangle */
1146  for (i = 0; i < xN; i++)
1147  for (j = i; j < xN; j++)
1148  y[j*yLD+i] = y[i*yLD+j];
1149 
1150  return SUCCESSFUL_RETURN;
1151 }
1152 
1153 
1155 
1156 
1157 /*
1158  * end of file
1159  */
int getMax(int x, int y)
Interfaces matrix-vector operations tailored to symmetric sparse matrices.
virtual returnValue transTimes(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
virtual BooleanType isDiag() const
virtual returnValue addToDiag(real_t alpha)
virtual returnValue getRow(int rNum, const Indexlist *const icols, real_t alpha, real_t *row) const
virtual Matrix * duplicate() const
virtual returnValue getCol(int cNum, const Indexlist *const irows, real_t alpha, real_t *col) const
BEGIN_NAMESPACE_ACADO const double EPS
virtual Matrix * duplicate() const
Allows to pass back messages to the calling function.
Interfaces matrix-vector operations tailored to symmetric dense matrices.
virtual returnValue getRow(int rNum, const Indexlist *const icols, real_t alpha, real_t *row) const
virtual returnValue getCol(int cNum, const Indexlist *const irows, real_t alpha, real_t *col) const
returnValue myPrintf(const char *s)
virtual Matrix * duplicate() const
virtual returnValue times(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
virtual real_t diag(int i) const
Interfaces matrix-vector operations tailored to general dense matrices.
Interfaces matrix-vector operations tailored to general sparse matrices.
virtual returnValue print() const
returnValue transTimes(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
virtual Matrix * duplicate() const
virtual returnValue print() const
virtual returnValue addToDiag(real_t alpha)
Abstract base class for interfacing tailored matrix-vector operations.
virtual returnValue bilinear(const Indexlist *const icols, int xN, const real_t *x, int xLD, real_t *y, int yLD) const
const char *const NOTRANS
virtual real_t diag(int i) const
BEGIN_NAMESPACE_QPOASES const char *const TRANS
#define v
BooleanType needToFreeMemory() const
#define BT_TRUE
Definition: acado_types.hpp:47
bool isExactlyMinusOne(real_t a)
RowXpr row(Index i)
Definition: BlockMethods.h:725
BEGIN_NAMESPACE_QPOASES bool isExactlyOne(real_t a)
returnValue times(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
virtual returnValue bilinear(const Indexlist *const icols, int xN, const real_t *x, int xLD, real_t *y, int yLD) const
BEGIN_NAMESPACE_QPOASES returnValue print(const real_t *const v, int n)
#define BT_FALSE
Definition: acado_types.hpp:49
virtual BooleanType isDiag() const
ColXpr col(Index i)
Definition: BlockMethods.h:708
double real_t
Definition: AD_test.c:10
bool isExactlyZero(real_t a)


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