dims_clapack.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2008, AIST, the University of Tokyo and General Robotix Inc.
3  * All rights reserved. This program is made available under the terms of the
4  * Eclipse Public License v1.0 which accompanies this distribution, and is
5  * available at http://www.eclipse.org/legal/epl-v10.html
6  * Contributors:
7  * The University of Tokyo
8  */
16 #include <dims_common.h>
17 #ifdef __darwin__
18 typedef int integer;
19 typedef double doublereal;
20 #include <stdlib.h>
21 #else
22 #include <f2c.h>
23 #include <malloc.h>
24 #endif
25 #include <stdio.h>
26 
27 #ifdef USE_CLAPACK_INTERFACE
28 #ifdef __WIN32__
29 #undef USE_CLAPACK_INTERFACE
30 #include <blaswrap.h>
31 #else
32 #if defined(__darwin__) || defined(__linux__)
33 extern "C"{
34 #endif
35 #include <cblas.h>
36 #if defined(__darwin__) || defined(__linux__)
37 }
38 #endif
39 #endif
40 #endif
41 
42 /* CLAPACK native functions */
43 extern "C"
44 {
45 int dgesvx_(
46  char *fact, char *trans, integer *n,
47  integer *nrhs, doublereal *a, integer *lda,
48  doublereal *af, integer *ldaf, integer *ipiv,
49  char *equed, doublereal *r, doublereal *c,
50  doublereal *b, integer *ldb, doublereal *x,
51  integer *ldx, doublereal *rcond, doublereal *ferr,
52  doublereal *berr, doublereal *work, integer *iwork,
53  integer *info);
54 
55 int dgesvd_(
56  char* jobu, char* jobvt, integer* m, integer* n,
57  doublereal *a, integer *lda, doublereal *s,
58  doublereal *u, integer *ldu, doublereal *vt,
59  integer *ldvt, doublereal* work, integer *lwork,
60  integer *info);
61 
62 int dgelss_(
63  integer* m, integer* n, integer* nrhs,
64  doublereal* A, integer* lda,
65  doublereal* B, integer* ldb,
66  doublereal* S, doublereal* rcond, integer* rank,
67  doublereal* work, integer* lwork, integer* info);
68 
69 int dporfs_(
70  char* uplo, integer* n, integer* nrhs,
71  doublereal* A, integer* lda,
72  doublereal* AF, integer* ldaf,
73  doublereal* b, integer* ldb,
74  doublereal* x, integer* ldx,
75  doublereal* ferr, doublereal* berr,
76  doublereal* work, integer* iwork,
77  integer* info);
78 
79 int dpotrf_(char* uplo, integer* n, doublereal* A, integer* lda,
80  integer* info);
81 
82 int dpotrs_(char* uplo, integer* n, integer* nrhs,
83  doublereal* A, integer* lda,
84  doublereal* x, integer* ldb,
85  integer* info);
86 
87 int dposv_(
88  char* uplo, integer* n, integer* nrhs,
89  doublereal* A, integer* lda,
90  doublereal* b, integer* ldb,
91  integer* info);
92 
93 int dposvx_(
94  char* fact, char* uplo, integer* n, integer* nrhs,
95  doublereal* A, integer* lda, doublereal* AF, integer* ldaf,
96  char* equed, doublereal* S,
97  doublereal* B, integer* ldb,
98  doublereal* X, integer* ldx,
99  doublereal* rcond, doublereal* ferr, doublereal* berr,
100  doublereal* work, integer* iwork, integer* info);
101 
102 int dgesvd_(
103  char* jobu, char* jobvt, integer* m, integer* n,
104  doublereal* A, integer* lda,
105  doublereal* S, doublereal* U, integer* ldu,
106  doublereal* VT, integer* ldvt,
107  doublereal* work, integer* lwork, integer* info);
108 
109 int dsyev_(char* jobz, char* uplo, integer* m, doublereal* a, integer* n, doublereal* w, doublereal* work, integer* lwork, integer* info);
110 
111 int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *a,
112  integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
113  integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
114  integer *lwork, integer *info);
115 
116 int dgetrf_(integer *m, integer *n, doublereal *a,
117  integer *lda, integer *ipiv, integer *info);
118 
119 #ifndef USE_CLAPACK_INTERFACE
120 int f2c_dscal(
121  integer* n, doublereal* alpha, doublereal* x, integer* incx);
122 
123 int f2c_dcopy(
124  integer* n, doublereal* x, integer* incx, doublereal* y, integer* incy);
125 double f2c_ddot(
126  integer* n, doublereal* x, integer* incx, doublereal* y, integer* incy);
127 int f2c_dgemv(
128  char* trans, integer* m, integer* n,
129  doublereal* alpha, doublereal* A, integer* lda,
130  doublereal* x, integer* incx,
131  doublereal* beta, doublereal* y, integer* incy);
132 
133 int f2c_dgemm(
134  char *transa, char *transb, integer *m, integer *n, integer *k,
135  doublereal *alpha, doublereal *a, integer *lda,
136  doublereal *b, integer *ldb,
137  doublereal *beta, doublereal *c, integer *ldc);
138 
139 int f2c_dsyrk(
140  char* uplo, char* trans, integer* n, integer* k,
141  doublereal* alpha, doublereal* A, integer* lda,
142  doublereal* beta, doublereal* C, integer* ldc);
143 
144 int f2c_daxpy(
145  integer* n, doublereal* alpha, doublereal* x, integer* incx,
146  doublereal* y, integer* incy);
147 #endif
148 }
149 
150 double dims_dot(double* _x, double* _y, int _n)
151 {
152 #ifndef USE_CLAPACK_INTERFACE
153  integer incx = 1, incy = 1;
154  integer n = _n;
155  return f2c_ddot(&n, _x, &incx, _y, &incy);
156 #else
157  const int n=_n, incX=1, incY=1;
158  return cblas_ddot(n, _x, incX, _y, incY);
159 #endif
160 
161 }
162 
163 int dims_copy(double* _x, double* _y, int _n)
164 {
165 #ifndef USE_CLAPACK_INTERFACE
166  integer incx = 1, incy = 1;
167  integer n = _n;
168  f2c_dcopy(&n, _x, &incx, _y, &incy);
169 #else
170  const int n = _n, incX = 1, incY = 1;
171  cblas_dcopy(n, _x, incX, _y, incY);
172 #endif
173  return 0;
174 }
175 
176 int dims_scale_myself(double* _x, double _alpha, int _n)
177 {
178 #ifndef USE_CLAPACK_INTERFACE
179  integer incx = 1;
180  integer n = _n;
181  f2c_dscal(&n, &_alpha, _x, &incx);
182 #else
183  const int n = _n, incX = 1;
184  const double alpha = _alpha;
185  cblas_dscal(n, alpha, _x, incX);
186 #endif
187  return 0;
188 }
189 
190 int dims_scale(double* _x, double _alpha, int _n, double*_y)
191 {
192 #ifndef USE_CLAPACK_INTERFACE
193  integer incx = 1, incy = 1;
194  integer n = _n;
195  f2c_dcopy(&n, _x, &incx, _y, &incy);
196  f2c_dscal(&n, &_alpha, _y, &incx);
197 #else
198  const int incX = 1, incY = 1;
199  cblas_dcopy(_n, _x, incX, _y, incY);
200  cblas_dscal(_n, _alpha, _x, incX);
201 #endif
202  return 0;
203 }
204 
205 int dims_dgemv(double* _A, int _m, int _n, double* _x, double* _y)
206 {
207 #ifndef USE_CLAPACK_INTERFACE
208  char trans = 'N';
209  integer m = _m;
210  integer n = _n;
211  integer lda = _m;
212  doublereal alpha = 1.0;
213  integer incx = 1;
214  doublereal beta = 0.0;
215  integer incy = 1;
216  f2c_dgemv(&trans, &m, &n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
217 #else
218  cblas_dgemv(CblasColMajor, CblasNoTrans,
219  _m, _n, 1.0, _A,_m, _x, 1, 0.0, _y, 1);
220 #endif
221  return 0;
222 }
223 
224 int dims_dgemv_tran(double* _A, int _m, int _n, double* _x, double* _y)
225 {
226 #ifndef USE_CLAPACK_INTERFACE
227  char trans = 'T';
228  integer m = _m;
229  integer n = _n;
230  integer lda = _m;
231  doublereal alpha = 1.0;
232  integer incx = 1;
233  doublereal beta = 0.0;
234  integer incy = 1;
235  f2c_dgemv(&trans, &m, &n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
236 #else
237  cblas_dgemv(CblasColMajor, CblasTrans,
238  _m, _n, 1.0, _A, _m, _x, 1, 0.0, _y, 1);
239 #endif
240  return 0;
241 }
242 
243 int dims_dgemm(double* _A, double* _B, int _m, int _n, int _k, double* _C)
244 {
245 #ifndef USE_CLAPACK_INTERFACE
246  char transa = 'N';
247  char transb = 'N';
248  integer m = _m;
249  integer n = _n;
250  integer k = _k;
251  doublereal alpha = 1.0;
252  integer lda = _m;
253  integer ldb = _k;
254  doublereal beta = 0.0;
255  integer ldc = _m;
256  f2c_dgemm(&transa, &transb, &m, &n, &k, &alpha, _A, &lda, _B, &ldb, &beta, _C, &ldc);
257 #else
258  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
259  _m, _n, _k, 1.0,
260  _A, _m, // lda
261  _B, _k, // ldb
262  0.0,
263  _C, _m // ldc
264  );
265 #endif
266  return 0;
267 }
268 
269 int dims_dsyrk(double* _A, int _n, int _k, double* _C)
270 {
271 #ifndef USE_CLAPACK_INTERFACE
272  char uplo = 'U';
273  char trans = 'N';
274  integer n = _n;
275  integer k = _k;
276  doublereal alpha = 1.0;
277  integer lda = _n;
278  doublereal beta = 0.0;
279  integer ldc = _n;
280 
281  f2c_dsyrk(&uplo, &trans, &n, &k, &alpha, _A, &lda, &beta, _C, &ldc);
282 #else
283  cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans,
284  _n, _k, 1.0, _A, _n, 0.0, _C, _n);
285 #endif
286  return 0;
287 }
288 
289 int dims_dsyrk_trans_first(double* _A, int _n, int _k, double* _C)
290 {
291 #ifndef USE_CLAPACK_INTERFACE
292  char uplo = 'U';
293  char trans = 'T';
294  integer n = _n;
295  integer k = _k;
296  doublereal alpha = 1.0;
297  integer lda = _k;
298  doublereal beta = 0.0;
299  integer ldc = _n;
300  f2c_dsyrk(&uplo, &trans, &n, &k, &alpha, _A, &lda, &beta, _C, &ldc);
301 #else
302  cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans,
303  _n, _k, 1.0, _A, _n, 0.0, _C, _n);
304 #endif
305  return 0;
306 }
307 
308 int dims_daxpy(int _n, double _alpha, double* _x, double* _y)
309 {
310 #ifndef USE_CLAPACK_INTERFACE
311  integer n = _n;
312  doublereal alpha = _alpha;
313  integer incx = 1;
314  integer incy = 1;
315  f2c_daxpy(&n, &alpha, _x, &incx, _y, &incy);
316 #else
317  cblas_daxpy(_n, _alpha, _x, 1, _y, 1);
318 #endif
319  return 0;
320 }
321 
322 int dims_dporfs(double* _a, double* _x, double* _b, int _m, int _nrhs)
323 {
324  char uplo = 'U';
325  integer n = _m;
326  integer nrhs = _nrhs;
327  integer lda = _m;
328  integer ldaf = _m;
329  integer ldb = _m;
330  integer ldx = _m;
331  integer info;
332  integer a_ndata = n*n;
333  integer b_ndata = n*nrhs;
334  integer i;
335  doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
336  doublereal* u = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
337  doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
338  doublereal* x = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
339  doublereal* ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
340  doublereal* berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
341  doublereal* work = (doublereal *)malloc(sizeof(doublereal)*3*n);
342  integer* iwork = (integer*)malloc(sizeof(integer)*n);
343  /* copy data in _a because they are overwritten by dposv() */
344  for(i=0; i<a_ndata; i++)
345  {
346  a[i] = _a[i];
347  u[i] = _a[i];
348  }
349  for(i=0; i<b_ndata; i++)
350  {
351  b[i] = _b[i];
352  x[i] = _b[i];
353  }
354  dpotrf_(&uplo, &n, u, &lda, &info);
355  dpotrs_(&uplo, &n, &nrhs, u, &lda, x, &ldb, &info);
356 #if 0
357  printf("info=%d\n", info);
358  printf("old x = ");
359  for(i=0; i<b_ndata; i++)
360  {
361  printf("\t%.8g", x[i]);
362  }
363  printf("\n");
364 #endif
365 /* dposv_(&uplo, &n, &nrhs, a, &lda, x, &ldb, &info);*/
366  dporfs_(&uplo, &n, &nrhs, a, &lda, u, &ldaf, b, &ldb, x, &ldx, ferr, berr, work, iwork, &info);
367  for(i=0; i<b_ndata; i++)
368  {
369  _x[i] = x[i];
370  }
371 #if 0
372  printf("info=%d\n", info);
373  printf("new x = ");
374  for(i=0; i<b_ndata; i++)
375  {
376  printf("\t%.8g", x[i]);
377  }
378  printf("\n");
379 #endif
380  free(a);
381  free(u);
382  free(b);
383  free(x);
384  free(ferr);
385  free(berr);
386  free(work);
387  free(iwork);
388  return info;
389 }
390 
391 int dims_dposvx(double* _a, double* _x, double* _b, int _m, int _nrhs, double* _rcond)
392 {
393  char fact = 'E';
394  char uplo = 'U';
395  integer n = _m;
396  integer nrhs = _nrhs;
397  integer lda = _m;
398  integer ldaf = _m;
399  char equed = 'N';
400  integer ldb = _m;
401  integer ldx = _m;
402  doublereal rcond;
403  integer info;
404  integer a_ndata = n*n;
405  integer b_ndata = n*nrhs;
406  integer i;
407  doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
408  doublereal* af = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
409  doublereal* s = (doublereal *)malloc(sizeof(doublereal)*n);
410  doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
411  doublereal* x = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
412  doublereal* ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
413  doublereal* berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
414  doublereal* work = (doublereal *)malloc(sizeof(doublereal)*3*n);
415  integer* iwork = (integer*)malloc(sizeof(integer)*n);
416  /* copy data in _a because they are overwritten by dposv() */
417  for(i=0; i<a_ndata; i++)
418  {
419  a[i] = _a[i];
420  }
421  for(i=0; i<b_ndata; i++)
422  {
423  b[i] = _b[i];
424  }
425  dposvx_(&fact, &uplo, &n, &nrhs, a, &lda, af, &ldaf, &equed, s, b, &ldb, x, &ldx, &rcond, ferr, berr, work, iwork, &info);
426  for(i=0; i<b_ndata; i++)
427  {
428  _x[i] = x[i];
429  }
430 #if 0
431  printf("scale=[");
432  for(i=0; i<n; i++)
433  {
434  printf("\t%.8g", s[i]);
435  }
436  printf("]\n");
437 #endif
438  *_rcond = rcond;
439  free(a);
440  free(af);
441  free(s);
442  free(b);
443  free(x);
444  free(ferr);
445  free(berr);
446  free(work);
447  free(iwork);
448  return info;
449 }
450 
451 int dims_dposv(double* _a, double* _x, double* _b, int _m, int _nrhs)
452 {
453  char uplo = 'U';
454  integer n = _m;
455  integer nrhs = _nrhs;
456  integer lda = _m;
457  integer ldb = _m;
458  integer info;
459  integer a_ndata = n*n;
460  integer b_ndata = n*nrhs;
461  integer i;
462  doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
463  doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
464  /* copy data in _a because they are overwritten by dposv() */
465  for(i=0; i<a_ndata; i++)
466  {
467  a[i] = _a[i];
468  }
469  for(i=0; i<b_ndata; i++)
470  {
471  b[i] = _b[i];
472  }
473  dposv_(&uplo, &n, &nrhs, a, &lda, b, &ldb, &info);
474  for(i=0; i<b_ndata; i++)
475  {
476  _x[i] = b[i];
477  }
478  free(a);
479  free(b);
480  return info;
481 }
482 
483 int dims_dgesvx(double* _a, double* _x, double* _b, int _n, int _nrhs)
484 {
485  char fact = 'N';
486  char trans = 'N';
487  integer n = (integer)_n;
488  integer nrhs = (integer)_nrhs;
489  doublereal* a = (doublereal*)_a;
490  integer lda = n;
491  doublereal* af = (doublereal*)malloc(sizeof(doublereal)*n*n);
492  integer ldaf = n;
493  integer *ipiv=(integer *)malloc(sizeof(integer)*n);
494  char equed = 'N';
495  doublereal *r=(doublereal *)malloc(sizeof(doublereal)*n);
496  doublereal *c=(doublereal *)malloc(sizeof(doublereal)*n);
497  doublereal* b = (doublereal*)_b;
498  integer ldb = n;
499  doublereal *x = (doublereal *)_x;
500  integer ldx = n;
501 
502  doublereal rcond;
503  doublereal *ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
504  doublereal *berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
505  doublereal *work = (doublereal *)malloc(sizeof(doublereal)*4*n);
506  integer *iwork = (integer *)malloc(sizeof(integer)*n);
507 
508  integer info;
509  dgesvx_(&fact, &trans, &n, &nrhs, a, &lda, af, &ldaf, ipiv,
510  &equed, r, c, b, &ldb, x, &ldx, &rcond, ferr, berr,
511  work, iwork, &info);
512  free(ipiv);
513  free(r);
514  free(c);
515  free(iwork);
516  free(work);
517  free(berr);
518  free(ferr);
519  free(af);
520  return info;
521 }
522 
523 /* inputs: A(m x n), b(m x nrhs), x(n x nrhs) */
524 int dims_dgelss(double* _a, double* _x, double* _b, int _m, int _n, int _nrhs,
525  double* _s, int* _rank, int _lwork)
526 {
527  int i, mm, tmp1, tmp2;
528  integer m = (integer)_m; /* row of A */
529  integer n = (integer)_n; /* col of A */
530  integer nrhs = (integer)_nrhs; /* row of b and x */
531  integer lda, ldb, rank, lwork, info;
532  doublereal *a, *b, *s, rcond, *work;
533  mm = _m*_n;
534  a = (doublereal*)malloc(sizeof(doublereal)*mm); /* A is overwritten by its singular vectors */
535  for(i=0; i<mm; i++) a[i] = _a[i];
536  lda = MAX(1, _m); /* row of A */
537  ldb = MAX3(1, _m, _n); /* row of B */
538  mm = ldb*_nrhs;
539  b = (doublereal*)malloc(sizeof(doublereal)*mm); /* note that b is overwritten by lapack function */
540  for(i=0; i<mm; i++) b[i] = _b[i]; /* copy b */
541  s = (doublereal*)_s;
542  rcond = -1; /* rcond < 0 -> use machine precision */
543  if(_lwork > 0) lwork = _lwork;
544  else
545  {
546  lwork = MIN(_m, _n); /* ? size of work */
547  lwork *= 3;
548  tmp1 = MIN(_m, _n);
549  tmp2 = MAX(_m, _n);
550  tmp1 *= 2;
551  lwork += MAX3(tmp1, tmp2, _nrhs);
552  lwork *= 5;
553  }
554  work = (doublereal*)malloc(sizeof(doublereal)*lwork);
555  dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank,
556  work, &lwork, &info);
557 
558  mm = _n*_nrhs;
559  for(i=0; i<mm; i++) _x[i] = b[i];
560  *_rank = (int)rank;
561 
562  free(a);
563  free(b);
564  free(work);
565  return info;
566 }
567 
568 int dims_svd(double* _a, int _m, int _n, double* _u, double* _sigma, double* _vt)
569 {
570  char jobu = 'A'; /* compute all columns of U */
571  char jobvt = 'A'; /* compute all rows of VT */
572  int mn = _m * _n, i, large = MAX(_m, _n), small = MIN(_m, _n);
573  integer m = (integer)_m;
574  integer n = (integer)_n;
575  integer lda = m;
576  integer ldu = m;
577  integer ldvt = n;
578  integer lwork, info;
579  doublereal* a, *work;
580  doublereal* u = (doublereal*)_u;
581  doublereal* s = (doublereal*)_sigma;
582  doublereal* vt = (doublereal*)_vt;
583  /* copy a - contents of a will be destroyed */
584  a = (doublereal*)malloc(sizeof(doublereal) * mn);
585  for(i=0; i<mn; i++) a[i] = _a[i];
586  /* compute lwork */
587  lwork = MAX(3*small+large, 5*small-4);
588  lwork *= 5; /* larger is better */
589  /* create work */
590  work = (doublereal*)malloc(sizeof(doublereal) * lwork);
591  /* call dgesvd */
592  dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);
593 
594  free(a);
595  free(work);
596  return info;
597 }
598 
599 
600 int dims_eigs(int _n, double *_a, double *w)
601 {
602  char jobz='V'; /* compute both eigenvalues and vectors */
603  char uplo='U'; /* store the upper triangle */
604  integer lwork;
605  double *work;
606  double *a;
607  integer info;
608  int i;
609  integer n = (integer)_n;
610 
611  a=(double*)malloc(sizeof(double)*n*n);
612  lwork=3*n;
613  lwork=lwork*5;
614  work=(double*)malloc(sizeof(double)*lwork);
615 
616  for(i=0;i<n*n;i++){
617  a[i]=_a[i];
618  }
619  dsyev_(&jobz,&uplo,&n,a,&n,w,work,&lwork,&info);
620  for(i=0;i<n*n;i++)
621  {
622  _a[i]=a[i];
623  }
624  free(a);
625  free(work);
626  return info;
627 }
628 
629 int dims_eigs2(int _n, double *_a, double *w)
630 {
631  char jobz='V'; /* computes both eigenvalues and vectors */
632  char uplo='U'; /* store the upper triangle */
633  integer lwork;
634  double *work=NULL;
635  double *a=NULL;
636  integer info;
637  int i;
638  integer n = (integer)_n;
639 
640  if(!a)
641  a=(double*)malloc(sizeof(double)*n*n);
642  lwork=3*n;
643  lwork=lwork*5;
644  if(!work)
645  work=(double*)malloc(sizeof(double)*lwork);
646 
647  for(i=0;i<n*n;i++)
648  {
649  a[i]=_a[i];
650  }
651  dsyev_(&jobz,&uplo,&n,a,&n,w,work,&lwork,&info);
652  for(i=0;i<n*n;i++)
653  {
654  _a[i]=a[i];
655  }
656  free(a);
657  free(work);
658 
659  return info;
660 }
661 
662 int dims_dgeev(int _n, double* _a, double* _wr, double* _wi, double* _vr)
663 {
664  int i, nn;
665  char jobvl = 'N';
666  char jobvr = 'V';
667  integer n, lda, ldvl, ldvr, lwork, info;
668  doublereal *a, *wr, *wi, *vl, *vr, *work;
669 
670  n = (integer)_n;
671  lda = n;
672  ldvl = 1;
673  ldvr = n;
674  lwork = 4*n;
675  nn = _n*_n;
676 
677  a = (doublereal*)malloc(sizeof(doublereal)*nn);
678  wr = (doublereal*)malloc(sizeof(doublereal)*n);
679  wi = (doublereal*)malloc(sizeof(doublereal)*n);
680  vl = (doublereal*)malloc(sizeof(doublereal)*ldvl*n);
681  vr = (doublereal*)malloc(sizeof(doublereal)*ldvr*n);
682  work = (doublereal *)malloc(sizeof(doublereal)*lwork);
683 
684  for(i=0; i<nn; i++) a[i] = _a[i];
685 
686  dgeev_(&jobvl, &jobvr, &n,
687  a, &lda,
688  wr, wi,
689  vl, &ldvl,
690  vr, &ldvr,
691  work, &lwork, &info);
692 
693  for(i=0; i<n; i++) _wr[i] = wr[i];
694  for(i=0; i<n; i++) _wi[i] = wi[i];
695  for(i=0; i<nn; i++) _vr[i] = vr[i];
696  free(a);
697  free(wr);
698  free(wi);
699  free(vl);
700  free(vr);
701  free(work);
702  return info;
703 }
704 
705 int dims_dgeev_simple(int _n, double* _a, double* _wr, double* _wi)
706 {
707  int i, nn;
708  char jobvl = 'N';
709  char jobvr = 'N';
710  integer n, lda, ldvl, ldvr, lwork, info;
711  doublereal *a, *wr, *wi, *vl, *vr, *work;
712 
713  n = (integer)_n;
714  lda = n;
715  ldvl = 1;
716  ldvr = n;
717  lwork = 3*n;
718  nn = _n*_n;
719 
720  a = (doublereal*)malloc(sizeof(doublereal)*nn);
721  wr = (doublereal*)malloc(sizeof(doublereal)*n);
722  wi = (doublereal*)malloc(sizeof(doublereal)*n);
723  vl = (doublereal*)malloc(sizeof(doublereal)*ldvl*n);
724  vr = (doublereal*)malloc(sizeof(doublereal)*ldvr*n);
725  work = (doublereal *)malloc(sizeof(doublereal)*lwork);
726 
727  for(i=0; i<nn; i++) a[i] = _a[i];
728 
729  dgeev_(&jobvl, &jobvr, &n,
730  a, &lda,
731  wr, wi,
732  vl, &ldvl,
733  vr, &ldvr,
734  work, &lwork, &info);
735 
736  for(i=0; i<n; i++) _wr[i] = wr[i];
737  for(i=0; i<n; i++) _wi[i] = wi[i];
738 
739  free(a);
740  free(wr);
741  free(wi);
742  free(vl);
743  free(vr);
744  free(work);
745  return info;
746 }
747 
748 int dims_det(int _n, double* _a, double* _x)
749 {
750  int i, nn, cnt = 0;
751  const double sign[2] = {1.0, -1.0};
752 
753  integer n, info;
754  doublereal *a;
755  integer* ipiv;
756 
757  n = (integer)_n;
758  nn = _n*_n;
759 
760  a = (doublereal*)malloc(sizeof(doublereal)*nn);
761  ipiv = (integer*)malloc(sizeof(integer)*n);
762 
763  for(i=0; i<nn; i++) a[i] = _a[i];
764 
765  dgetrf_(&n, &n,
766  a, &n,
767  ipiv, &info);
768 
769  *_x = 1.0;
770  for(i=0; i<n; i++)
771  {
772  *_x *= a[i*n+i];
773  cnt += ((i+1)!=ipiv[i]);
774  }
775  *_x *= sign[cnt%2];
776 /* if(cnt%2 != 0) *_x *= -1.0; */
777  free(a);
778  free(ipiv);
779  return info;
780 }
int dims_dgeev_simple(int _n, double *_a, double *_wr, double *_wi)
Computes eigenvalues only.
int c
Definition: autoplay.py:16
#define MAX(a, b)
Returns the max value between a and b.
Definition: IceTypes.h:147
int dims_dgemv(double *_A, int _m, int _n, double *_x, double *_y)
* x
Definition: IceUtils.h:98
int dims_scale(double *_x, double _alpha, int _n, double *_y)
int dims_dgemv_tran(double *_A, int _m, int _n, double *_x, double *_y)
Definition: clap.cpp:13
int dsyev_(char *jobz, char *uplo, integer *m, doublereal *a, integer *n, doublereal *w, doublereal *work, integer *lwork, integer *info)
int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info)
double dims_dot(double *_x, double *_y, int _n)
char * malloc()
int dims_det(int _n, double *_a, double *_x)
Computes the determinant.
#define MIN(a, b)
Returns the min value between a and b.
Definition: IceTypes.h:146
int dims_dgeev(int _n, double *_a, double *_wr, double *_wi, double *_vr)
Computes eigenvalues and eigenvectors.
int dims_eigs(int _n, double *_a, double *w)
Eigenvalues / eigenvectors.
png_uint_32 i
Definition: png.h:2735
long b
Definition: jpegint.h:371
png_infop png_bytep * trans
Definition: png.h:2435
int dims_dposvx(double *_a, double *_x, double *_b, int _m, int _nrhs, double *_rcond)
int dpotrf_(char *uplo, integer *n, doublereal *A, integer *lda, integer *info)
int dporfs_(char *uplo, integer *n, integer *nrhs, doublereal *A, integer *lda, doublereal *AF, integer *ldaf, doublereal *b, integer *ldb, doublereal *x, integer *ldx, doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork, integer *info)
#define MAX3(l, m, n)
Definition: dims_common.h:29
int dpotrs_(char *uplo, integer *n, integer *nrhs, doublereal *A, integer *lda, doublereal *x, integer *ldb, integer *info)
string a
int dims_eigs2(int _n, double *_a, double *w)
int free()
int dims_daxpy(int _n, double _alpha, double *_x, double *_y)
int f2c_daxpy(integer *n, doublereal *alpha, doublereal *x, integer *incx, doublereal *y, integer *incy)
int f2c_dcopy(integer *n, doublereal *x, integer *incx, doublereal *y, integer *incy)
int dims_dporfs(double *_a, double *_x, double *_b, int _m, int _nrhs)
For positive-definite, symmetric matrices.
int dims_dgesvx(double *_a, double *_x, double *_b, int _n, int _nrhs)
Solves linear equation using LU decomposition.
int dims_dgemm(double *_A, double *_B, int _m, int _n, int _k, double *_C)
int dposvx_(char *fact, char *uplo, integer *n, integer *nrhs, doublereal *A, integer *lda, doublereal *AF, integer *ldaf, char *equed, doublereal *S, doublereal *B, integer *ldb, doublereal *X, integer *ldx, doublereal *rcond, doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork, integer *info)
int f2c_dscal(integer *n, doublereal *alpha, doublereal *x, integer *incx)
typedef int
Definition: png.h:1113
backing_store_ptr info
Definition: jmemsys.h:181
int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n, doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, integer *info)
int dims_dsyrk_trans_first(double *_A, int _n, int _k, double *_C)
int f2c_dgemv(char *trans, integer *m, integer *n, doublereal *alpha, doublereal *A, integer *lda, doublereal *x, integer *incx, doublereal *beta, doublereal *y, integer *incy)
double f2c_ddot(integer *n, doublereal *x, integer *incx, doublereal *y, integer *incy)
int dims_scale_myself(double *_x, double _alpha, int _n)
int dims_dposv(double *_a, double *_x, double *_b, int _m, int _nrhs)
* y
Definition: IceUtils.h:97
int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work, integer *lwork, integer *info)
int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *A, integer *lda, doublereal *B, integer *ldb, doublereal *S, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, integer *info)
int dims_dsyrk(double *_A, int _n, int _k, double *_C)
int dims_copy(double *_x, double *_y, int _n)
Wrappers of BLAS functions.
int dposv_(char *uplo, integer *n, integer *nrhs, doublereal *A, integer *lda, doublereal *b, integer *ldb, integer *info)
int f2c_dgemm(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
int dims_dgelss(double *_a, double *_x, double *_b, int _m, int _n, int _nrhs, double *_s, int *_rank, int _lwork)
Solves linear equation using singular-value decomposition.
int dgesvx_(char *fact, char *trans, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *af, integer *ldaf, integer *ipiv, char *equed, doublereal *r, doublereal *c, doublereal *b, integer *ldb, doublereal *x, integer *ldx, doublereal *rcond, doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork, integer *info)
int dims_svd(double *_a, int _m, int _n, double *_u, double *_sigma, double *_vt)
Performs singular value decomposition.
int f2c_dsyrk(char *uplo, char *trans, integer *n, integer *k, doublereal *alpha, doublereal *A, integer *lda, doublereal *beta, doublereal *C, integer *ldc)


openhrp3
Author(s): AIST, General Robotix Inc., Nakamura Lab of Dept. of Mechano Informatics at University of Tokyo
autogenerated on Sat Apr 13 2019 02:14:21