40 for (i = 0; i <
N; i++) av[i] = i+1;
41 for (i = 0; i <
N; i++) aTv[i] = i+1;
42 for (i = 0; i <
N; i++) bv[i] = i+1;
49 a.times(1, 1.0, bv, N, 0.0, &c, 1);
50 fprintf(
stdFile,
"Dot product; Error in sum of first %d squares: %9.2e\n", N, c - (1.0/6.0)*N*(N+1)*(2*N+1));
51 aT.transTimes(1, 1.0, bv, N, 0.0, &c, 1);
52 fprintf(
stdFile,
"Transpose; Error in sum of first %d squares: %9.2e\n", N, c - (1.0/6.0)*N*(N+1)*(2*N+1));
63 real_t _Av[] = {1.0/3.0, 1.0, 0.5, 0.25,
64 0.25, 0.5, 1.0/3.0, 0.2,
65 0.2, 1.0/3.0, 0.25, 1.0/6.0,
66 1.0/6.0, 0.25, 0.2, 1.0/7.0};
68 real_t Bv[] = {240, 16, -120, -140,
69 -2700, -120, 1200, 1680,
70 6480, 240, -2700, -4200,
71 -4200, -140, 1680, 2800};
80 for (i = 0; i < 16; i++) Av[i] = _Av[i];
82 A.times(4, 1.0, Bv, 4, 0.0, Cv, 4);
85 for (j = 0; j < 4; j++)
87 for (i = 0; i < 4; i++)
89 d = fabs(Cv[j*4+i] - static_cast<real_t>(i == j));
93 fprintf(
stdFile,
"Hilbert; Deviation from identity: %9.2e\n", err);
107 real_t _Asubv[] = {1.0/3.0, 0.25,
111 real_t Bsubv[] = {240, 16, -120, -140,
112 -2700, -120, 1200, 1680};
117 for (i = 0; i < 8; i++) Asubv[i] = _Asubv[i];
122 Asub.transTimes(2, 1.0, Bsubv, 4, 0.0, Csubv, 2);
125 for (j = 0; j < 2; j++)
127 for (i = 0; i < 2; i++)
129 d = fabs(Csubv[j*2+i] - static_cast<real_t>(i == j));
130 if (d > err) err = d;
133 fprintf(
stdFile,
"Submatrix transpose; Deviation from identity: %9.2e\n", err);
143 const int M = 20,
N = 15, K = 5;
148 real_t *Av, *X, *Y, *x, *
y, *xc, *yc;
153 for (i = 0; i < m; i++) { rows.addNumber(M-1 - 2*i); rows.addNumber(2*i); }
158 rows.getNumberArray(&rNum);
160 for (i = 0; i < m; i++) fprintf(
stdFile,
" %2d", rNum[i]);
165 for (i = 0; i < n; i++) fprintf(
stdFile,
" %2d", cNum[i]);
179 for (i = 0; i < M*
N; i++) Av[i] = -0.5*N*M + (
real_t)i;
180 for (i = 0; i < n*K; i++) X[i] = 1.0 / (
real_t)(i+1);
181 for (i = 0; i < m*K; i++) Y[i] = 1.0 / (
real_t)(i+1);
184 A.times(&rows, &cols, K, 1.0, X, n, 0.0, y, m);
187 for (j = 0; j < m; j++)
189 for (k = 0; k < K; k++)
191 yc[j+k*m] = -y[j+k*m];
192 for (i = 0; i < n; i++)
193 yc[j+k*m] += Av[cNum[i]+rNum[j]*N] * X[i+k*n];
194 if (fabs(yc[j+k*m]) > err) err = fabs(yc[j+k*m]);
197 fprintf(
stdFile,
"Indexlist submatrix; error: %9.2e\n", err);
200 A.transTimes(&rows, &cols, K, 1.0, Y, m, 0.0, x, n);
204 for (j = 0; j < n; j++)
206 for (k = 0; k < K; k++)
208 xc[j+k*n] = -x[j+k*n];
209 for (i = 0; i < m; i++)
210 xc[j+k*n] += Av[cNum[j]+rNum[i]*N] * Y[i+k*m];
211 if (fabs(xc[j+k*n]) > err) err = fabs(xc[j+k*n]);
214 fprintf(
stdFile,
"Indexlist transpose submatrix; error: %9.2e\n", err);
247 jc[0] = 0; jc[1] = 2; jc[2] = 5; jc[3] = 10;
248 ir[0] = 1; ir[1] = 2;
249 ir[2] = 0; ir[3] = 3; ir[4] = 4;
250 ir[5] = 0; ir[6] = 1; ir[7] = 2; ir[8] = 3; ir[9] = 4;
251 for (i = 0; i < 10; i++) val[i] = 1.0+i;
270 for (j = 0; j < 3; j++)
272 fprintf(
stdFile,
"Column %ld:\n", j);
273 A.
getCol(j, &rows, 1.0, col);
274 for (i = 0; i < 4; i++)
275 fprintf(
stdFile,
" %3.0f\n", col[i]);
299 jc[0] = 0; jc[1] = 2; jc[2] = 3; jc[3] = 5; jc[4] = 8; jc[5] = 10;
300 ir[0] = 1; ir[1] = 2;
302 ir[3] = 0; ir[4] = 2;
303 ir[5] = 0; ir[6] = 1; ir[7] = 2;
304 ir[8] = 1; ir[9] = 2;
305 for (i = 0; i < 10; i++) val[i] = 1.0+i;
323 for (j = 0; j < 3; j++)
325 A.
getRow(j, &cols, 1.0, row);
326 for (i = 0; i < 4; i++)
327 fprintf(
stdFile,
" %3.0f", row[i]);
347 real_t Ax[] = {-23, -11, -26, 42, 74, 99};
348 real_t ATy[] = {-63, -69, -222, -423, -359, 272, 126, 663, 1562, 1656};
351 for (i = 0; i < 10; i++) x[i] = -4.0+i;
360 jc[0] = 0; jc[1] = 2; jc[2] = 3; jc[3] = 5; jc[4] = 8; jc[5] = 10;
361 ir[0] = 1; ir[1] = 2;
363 ir[3] = 0; ir[4] = 2;
364 ir[5] = 0; ir[6] = 1; ir[7] = 2;
365 ir[8] = 1; ir[9] = 2;
366 for (i = 0; i < 10; i++) val[i] = 1.0+i;
370 A.
times(2, 1.0, x, 5, 0.0, y, 3);
372 for (i = 0; i < 6; i++)
373 if (fabs(y[i] - Ax[i]) > err) err = fabs(y[i] - Ax[i]);
374 fprintf(
stdFile,
"Error in sparse A*x: %9.2e\n", err);
379 for (i = 0; i < 10; i++)
380 if (fabs(x[i] - ATy[i]) > err) err = fabs(x[i] - ATy[i]);
381 fprintf(
stdFile,
"Error in sparse A'*x: %9.2e\n", err);
394 long nRows = 2 * N + 1;
402 real_t Ax[] = {0.31, 0.05, 0.06, 0.30, 0.76, 0.20, 0.24, 0.60};
403 real_t ATy[] = {0.278, 0.000, 0.548, 0.776, 0.000, 1.208};
405 Indexlist rows(4), cols(3), allcols(nCols);
416 for (i = 0; i < nCols; i++)
420 for (i = 0; i <=
N; i++) jc[i] = 3*i;
421 for (j = 0; j <
N; j++)
422 for (i = 0; i < 3; i++)
425 val[j*3+i] = 1.0 - 0.1 * (j*3+i);
429 fprintf(
stdFile,
"Test matrix A =\n");
430 for (j = 0; j < nRows; j++)
432 A.
getRow(j, &allcols, 1.0, xc);
433 for (i = 0; i < nCols; i++)
434 fprintf(
stdFile,
"%6.2f", xc[i]);
438 for (i = 0; i < 6; i++)
439 xc[i] = (1.0 + i) * 0.1;
441 A.
times(&rows, &cols, 2, 1.0, xc, 3, 0.0, yc, 4,
BT_TRUE);
444 for (i = 0; i < 8; i++)
445 if (fabs(yc[i] - Ax[i]) > err)
446 err = fabs(yc[i] - Ax[i]);
447 fprintf(
stdFile,
"Error in sparse indexed A*x: %9.2e\n", err);
449 A.
transTimes(&rows, &cols, 2, 1.0, yc, 4, 0.0, xc, 3);
451 for (i = 0; i < 6; i++)
452 if (fabs(xc[i] - ATy[i]) > err)
453 err = fabs(xc[i] - ATy[i]);
454 fprintf(
stdFile,
"Error in sparse indexed A'*y: %9.2e\n", err);
480 jc[0] = 0; jc[1] = 2; jc[2] = 5; jc[3] = 10;
481 ir[0] = 1; ir[1] = 2;
482 ir[2] = 0; ir[3] = 3; ir[4] = 4;
483 ir[5] = 0; ir[6] = 1; ir[7] = 2; ir[8] = 3; ir[9] = 4;
484 for (i = 0; i < 10; i++) val[i] = 1.0+i;
507 for (j = 0; j < 3; j++)
509 fprintf(
stdFile,
"Column %ld:\n", j);
510 A.
getCol(j, &rows, 1.0, col);
511 for (i = 0; i < 4; i++)
512 fprintf(
stdFile,
" %3.0f\n", col[i]);
536 jc[0] = 0; jc[1] = 2; jc[2] = 3; jc[3] = 5; jc[4] = 8; jc[5] = 10;
537 ir[0] = 1; ir[1] = 2;
539 ir[3] = 0; ir[4] = 2;
540 ir[5] = 0; ir[6] = 1; ir[7] = 2;
541 ir[8] = 1; ir[9] = 2;
542 for (i = 0; i < 10; i++) val[i] = 1.0+i;
564 for (j = 0; j < 3; j++)
566 A.
getRow(j, &cols, 1.0, row);
567 for (i = 0; i < 4; i++)
568 fprintf(
stdFile,
" %3.0f", row[i]);
588 real_t Ax[] = {-23, -11, -26, 42, 74, 99};
589 real_t ATy[] = {-63, -69, -222, -423, -359, 272, 126, 663, 1562, 1656};
592 for (i = 0; i < 10; i++) x[i] = -4.0+i;
601 jc[0] = 0; jc[1] = 2; jc[2] = 3; jc[3] = 5; jc[4] = 8; jc[5] = 10;
602 ir[0] = 1; ir[1] = 2;
604 ir[3] = 0; ir[4] = 2;
605 ir[5] = 0; ir[6] = 1; ir[7] = 2;
606 ir[8] = 1; ir[9] = 2;
607 for (i = 0; i < 10; i++) val[i] = 1.0+i;
615 A.
times(2, 1.0, x, 5, 0.0, y, 3);
617 for (i = 0; i < 6; i++)
618 if (fabs(y[i] - Ax[i]) > err) err = fabs(y[i] - Ax[i]);
619 fprintf(
stdFile,
"Error in sparse A*x: %9.2e\n", err);
624 for (i = 0; i < 10; i++)
625 if (fabs(x[i] - ATy[i]) > err) err = fabs(x[i] - ATy[i]);
626 fprintf(
stdFile,
"Error in sparse A'*x: %9.2e\n", err);
639 long nRows = 2 * N + 1;
647 real_t Ax[] = {0.31, 0.05, 0.06, 0.30, 0.76, 0.20, 0.24, 0.60};
648 real_t ATy[] = {0.278, 0.000, 0.548, 0.776, 0.000, 1.208};
650 Indexlist rows(4), cols(3), allcols(nCols);
661 for (i = 0; i < nCols; i++)
665 for (i = 0; i <=
N; i++) jc[i] = 3*i;
666 for (j = 0; j <
N; j++)
667 for (i = 0; i < 3; i++)
670 val[j*3+i] = 1.0 - 0.1 * (j*3+i);
678 fprintf(
stdFile,
"Test matrix A =\n");
679 for (j = 0; j < nRows; j++)
681 A.
getRow(j, &allcols, 1.0, xc);
682 for (i = 0; i < nCols; i++)
683 fprintf(
stdFile,
"%6.2f", xc[i]);
687 for (i = 0; i < 6; i++)
688 xc[i] = (1.0 + i) * 0.1;
690 A.
times(&rows, &cols, 2, 1.0, xc, 3, 0.0, yc, 4,
BT_TRUE);
693 for (i = 0; i < 8; i++)
694 if (fabs(yc[i] - Ax[i]) > err)
695 err = fabs(yc[i] - Ax[i]);
696 fprintf(
stdFile,
"Error in sparse indexed A*x: %9.2e\n", err);
698 A.
transTimes(&rows, &cols, 2, 1.0, yc, 4, 0.0, xc, 3);
699 for (i = 0; i < 3; i++)
701 for (j = 0; j < 2; j++)
702 fprintf(
stdFile,
"%6.2f", ATy[i + j*3]);
705 for (i = 0; i < 3; i++)
707 for (j = 0; j < 2; j++)
708 fprintf(
stdFile,
"%6.2f", xc[i + j*3]);
712 for (i = 0; i < 6; i++)
713 if (fabs(xc[i] - ATy[i]) > err)
714 err = fabs(xc[i] - ATy[i]);
715 fprintf(
stdFile,
"Error in sparse indexed A'*y: %9.2e\n", err);
731 real_t ZHZv[] = {0.144, 0.426, 0.708, 0.426, 1.500, 2.574, 0.708, 2.574, 4.440};
737 for (i = 0; i < 36; i++) Hv[i] = 0.0;
738 for (i = 0; i < 6; i++) Hv[i*7] = 1.0 - 0.1 * i;
739 for (i = 0; i < 5; i++) Hv[i*7+1] = Hv[i*7+6] = -0.1 * (i+1);
745 for (i = 0; i < 6; ++i)
747 for (j = 0; j < 6; ++j)
748 fprintf (stderr,
"%3.3f ", Hv[i*6+j]);
749 fprintf (stderr,
"\n");
751 fprintf (stderr,
"\n");
757 for (i = 0; i < 18; i++) Z[i] = 0.1 * (i+1);
759 fprintf (stderr,
"\n");
760 for (i = 0; i < 6; ++i)
762 for (j = 0; j < 3; ++j)
763 fprintf (stderr,
"%3.3f ", Z[i+j*6]);
764 fprintf (stderr,
"\n");
766 fprintf (stderr,
"\n");
768 Hd->
bilinear(cols, 3, Z, 6, ZHZd, 3);
771 for (i = 0; i < 9; i++)
772 if (fabs(ZHZd[i] - ZHZv[i]) > err)
773 err = fabs(ZHZd[i] - ZHZv[i]);
774 fprintf(
stdFile,
"Error in indexed dense bilinear form: %9.2e\n", err);
776 Hs->
bilinear(cols, 3, Z, 6, ZHZs, 3);
778 for (i = 0; i < 3; ++i)
780 for (j = 0; j < 3; ++j)
781 fprintf (stderr,
"%3.3f ", ZHZd[i*3+j]);
782 fprintf (stderr,
"\n");
784 fprintf (stderr,
"\n");
786 for (i = 0; i < 3; ++i)
788 for (j = 0; j < 3; ++j)
789 fprintf (stderr,
"%3.3f ", ZHZv[i*3+j]);
790 fprintf (stderr,
"\n");
792 fprintf (stderr,
"\n");
795 for (i = 0; i < 9; i++)
796 if (fabs(ZHZs[i] - ZHZv[i]) > err)
797 err = fabs(ZHZs[i] - ZHZv[i]);
798 fprintf(
stdFile,
"Error in indexed sparse bilinear form: %9.2e\n", err);
USING_NAMESPACE_QPOASES int sumOfSquares()
#define USING_NAMESPACE_QPOASES
returnValue addNumber(int addnumber)
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
returnValue getNumberArray(int *const numberarray) const
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
int indexDenseSubmatrix()
virtual returnValue getCol(int cNum, const Indexlist *const irows, real_t alpha, real_t *col) const
virtual returnValue transTimes(int_t xN, real_t alpha, const real_t *x, int_t xLD, real_t beta, real_t *y, int_t yLD) const
Interfaces matrix-vector operations tailored to general sparse matrices.
virtual returnValue times(int xN, real_t alpha, const real_t *x, int xLD, real_t beta, real_t *y, int yLD) const
virtual returnValue getRow(int_t rNum, const Indexlist *const icols, real_t alpha, real_t *row) const
virtual returnValue times(int_t xN, real_t alpha, const real_t *x, int_t xLD, real_t beta, real_t *y, int_t yLD) const
void DenseMatrixCON(DenseMatrix *_THIS, int m, int n, int lD, real_t *v)
Interfaces matrix-vector operations tailored to general dense matrices.
Interfaces matrix-vector operations tailored to general sparse matrices.
virtual returnValue getCol(int_t cNum, const Indexlist *const irows, real_t alpha, real_t *col) const
virtual returnValue bilinear(const Indexlist *const icols, int xN, const real_t *x, int xLD, real_t *y, int yLD) const
Stores and manages index lists.
virtual returnValue bilinear(const Indexlist *const icols, int xN, const real_t *x, int xLD, real_t *y, int yLD) const