19 typedef double doublereal;
27 #ifdef USE_CLAPACK_INTERFACE
29 #undef USE_CLAPACK_INTERFACE
32 #if defined(__darwin__) || defined(__linux__)
36 #if defined(__darwin__) || defined(__linux__)
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,
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,
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);
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,
79 int dpotrf_(
char* uplo, integer*
n, doublereal* A, integer* lda,
82 int dpotrs_(
char* uplo, integer*
n, integer* nrhs,
83 doublereal* A, integer* lda,
84 doublereal* x, integer* ldb,
88 char* uplo, integer*
n, integer* nrhs,
89 doublereal* A, integer* lda,
90 doublereal*
b, integer* ldb,
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);
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);
109 int dsyev_(
char* jobz,
char* uplo, integer* m, doublereal*
a, integer*
n, doublereal* w, doublereal* work, integer* lwork, integer*
info);
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);
116 int dgetrf_(integer *m, integer *
n, doublereal *
a,
117 integer *lda, integer *ipiv, integer *
info);
119 #ifndef USE_CLAPACK_INTERFACE
121 integer*
n, doublereal* alpha, doublereal* x, integer* incx);
124 integer*
n, doublereal* x, integer* incx, doublereal* y, integer* incy);
126 integer*
n, doublereal* x, integer* incx, doublereal* y, integer* incy);
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);
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);
140 char* uplo,
char*
trans, integer*
n, integer* k,
141 doublereal* alpha, doublereal* A, integer* lda,
142 doublereal* beta, doublereal* C, integer* ldc);
145 integer*
n, doublereal* alpha, doublereal* x, integer* incx,
146 doublereal* y, integer* incy);
152 #ifndef USE_CLAPACK_INTERFACE
153 integer incx = 1, incy = 1;
155 return f2c_ddot(&
n, _x, &incx, _y, &incy);
157 const int n=_n, incX=1, incY=1;
158 return cblas_ddot(
n, _x, incX, _y, incY);
165 #ifndef USE_CLAPACK_INTERFACE
166 integer incx = 1, incy = 1;
170 const int n = _n, incX = 1, incY = 1;
171 cblas_dcopy(
n, _x, incX, _y, incY);
178 #ifndef USE_CLAPACK_INTERFACE
183 const int n = _n, incX = 1;
184 const double alpha = _alpha;
185 cblas_dscal(
n, alpha, _x, incX);
192 #ifndef USE_CLAPACK_INTERFACE
193 integer incx = 1, incy = 1;
198 const int incX = 1, incY = 1;
199 cblas_dcopy(_n, _x, incX, _y, incY);
200 cblas_dscal(_n, _alpha, _x, incX);
205 int dims_dgemv(
double* _A,
int _m,
int _n,
double* _x,
double* _y)
207 #ifndef USE_CLAPACK_INTERFACE
212 doublereal alpha = 1.0;
214 doublereal beta = 0.0;
216 f2c_dgemv(&
trans, &m, &
n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
218 cblas_dgemv(CblasColMajor, CblasNoTrans,
219 _m, _n, 1.0, _A,_m, _x, 1, 0.0, _y, 1);
226 #ifndef USE_CLAPACK_INTERFACE
231 doublereal alpha = 1.0;
233 doublereal beta = 0.0;
235 f2c_dgemv(&
trans, &m, &
n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
237 cblas_dgemv(CblasColMajor, CblasTrans,
238 _m, _n, 1.0, _A, _m, _x, 1, 0.0, _y, 1);
243 int dims_dgemm(
double* _A,
double* _B,
int _m,
int _n,
int _k,
double* _C)
245 #ifndef USE_CLAPACK_INTERFACE
251 doublereal alpha = 1.0;
254 doublereal beta = 0.0;
256 f2c_dgemm(&transa, &transb, &m, &
n, &k, &alpha, _A, &lda, _B, &ldb, &beta, _C, &ldc);
258 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
271 #ifndef USE_CLAPACK_INTERFACE
276 doublereal alpha = 1.0;
278 doublereal beta = 0.0;
283 cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans,
284 _n, _k, 1.0, _A, _n, 0.0, _C, _n);
291 #ifndef USE_CLAPACK_INTERFACE
296 doublereal alpha = 1.0;
298 doublereal beta = 0.0;
302 cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans,
303 _n, _k, 1.0, _A, _n, 0.0, _C, _n);
308 int dims_daxpy(
int _n,
double _alpha,
double* _x,
double* _y)
310 #ifndef USE_CLAPACK_INTERFACE
312 doublereal alpha = _alpha;
317 cblas_daxpy(_n, _alpha, _x, 1, _y, 1);
322 int dims_dporfs(
double* _a,
double* _x,
double* _b,
int _m,
int _nrhs)
326 integer nrhs = _nrhs;
332 integer a_ndata =
n*
n;
333 integer b_ndata =
n*nrhs;
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);
344 for(
i=0;
i<a_ndata;
i++)
349 for(
i=0;
i<b_ndata;
i++)
357 printf(
"info=%d\n",
info);
359 for(
i=0;
i<b_ndata;
i++)
361 printf(
"\t%.8g", x[
i]);
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++)
372 printf(
"info=%d\n",
info);
374 for(
i=0;
i<b_ndata;
i++)
376 printf(
"\t%.8g", x[
i]);
391 int dims_dposvx(
double* _a,
double* _x,
double* _b,
int _m,
int _nrhs,
double* _rcond)
396 integer nrhs = _nrhs;
404 integer a_ndata =
n*
n;
405 integer b_ndata =
n*nrhs;
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);
417 for(
i=0;
i<a_ndata;
i++)
421 for(
i=0;
i<b_ndata;
i++)
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++)
434 printf(
"\t%.8g", s[
i]);
451 int dims_dposv(
double* _a,
double* _x,
double* _b,
int _m,
int _nrhs)
455 integer nrhs = _nrhs;
459 integer a_ndata =
n*
n;
460 integer b_ndata =
n*nrhs;
462 doublereal*
a = (doublereal *)
malloc(
sizeof(doublereal)*a_ndata);
463 doublereal*
b = (doublereal *)
malloc(
sizeof(doublereal)*b_ndata);
465 for(
i=0;
i<a_ndata;
i++)
469 for(
i=0;
i<b_ndata;
i++)
474 for(
i=0;
i<b_ndata;
i++)
483 int dims_dgesvx(
double* _a,
double* _x,
double* _b,
int _n,
int _nrhs)
487 integer
n = (integer)_n;
488 integer nrhs = (integer)_nrhs;
489 doublereal*
a = (doublereal*)_a;
491 doublereal* af = (doublereal*)
malloc(
sizeof(doublereal)*
n*
n);
493 integer *ipiv=(integer *)
malloc(
sizeof(integer)*
n);
495 doublereal *r=(doublereal *)
malloc(
sizeof(doublereal)*
n);
496 doublereal *
c=(doublereal *)
malloc(
sizeof(doublereal)*
n);
497 doublereal*
b = (doublereal*)_b;
499 doublereal *x = (doublereal *)_x;
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);
510 &equed, r,
c,
b, &ldb, x, &ldx, &rcond, ferr, berr,
524 int dims_dgelss(
double* _a,
double* _x,
double* _b,
int _m,
int _n,
int _nrhs,
525 double* _s,
int* _rank,
int _lwork)
527 int i, mm, tmp1, tmp2;
528 integer m = (integer)_m;
529 integer
n = (integer)_n;
530 integer nrhs = (integer)_nrhs;
531 integer lda, ldb, rank, lwork,
info;
532 doublereal *
a, *
b, *s, rcond, *work;
534 a = (doublereal*)
malloc(
sizeof(doublereal)*mm);
535 for(
i=0;
i<mm;
i++)
a[
i] = _a[
i];
537 ldb =
MAX3(1, _m, _n);
539 b = (doublereal*)
malloc(
sizeof(doublereal)*mm);
540 for(
i=0;
i<mm;
i++)
b[
i] = _b[
i];
543 if(_lwork > 0) lwork = _lwork;
551 lwork +=
MAX3(tmp1, tmp2, _nrhs);
554 work = (doublereal*)
malloc(
sizeof(doublereal)*lwork);
555 dgelss_(&m, &
n, &nrhs,
a, &lda,
b, &ldb, s, &rcond, &rank,
556 work, &lwork, &
info);
559 for(
i=0;
i<mm;
i++) _x[
i] =
b[
i];
568 int dims_svd(
double* _a,
int _m,
int _n,
double* _u,
double* _sigma,
double* _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;
579 doublereal*
a, *work;
580 doublereal* u = (doublereal*)_u;
581 doublereal* s = (doublereal*)_sigma;
582 doublereal* vt = (doublereal*)_vt;
584 a = (doublereal*)
malloc(
sizeof(doublereal) * mn);
585 for(
i=0;
i<mn;
i++)
a[
i] = _a[
i];
587 lwork =
MAX(3*small+large, 5*small-4);
590 work = (doublereal*)
malloc(
sizeof(doublereal) * lwork);
592 dgesvd_(&jobu, &jobvt, &m, &
n,
a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &
info);
609 integer
n = (integer)_n;
614 work=(
double*)
malloc(
sizeof(
double)*lwork);
638 integer
n = (integer)_n;
645 work=(
double*)
malloc(
sizeof(
double)*lwork);
662 int dims_dgeev(
int _n,
double* _a,
double* _wr,
double* _wi,
double* _vr)
667 integer
n, lda, ldvl, ldvr, lwork,
info;
668 doublereal *
a, *wr, *wi, *vl, *vr, *work;
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);
684 for(
i=0;
i<nn;
i++)
a[
i] = _a[
i];
691 work, &lwork, &
info);
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];
710 integer
n, lda, ldvl, ldvr, lwork,
info;
711 doublereal *
a, *wr, *wi, *vl, *vr, *work;
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);
727 for(
i=0;
i<nn;
i++)
a[
i] = _a[
i];
734 work, &lwork, &
info);
736 for(
i=0;
i<
n;
i++) _wr[
i] = wr[
i];
737 for(
i=0;
i<
n;
i++) _wi[
i] = wi[
i];
751 const double sign[2] = {1.0, -1.0};
760 a = (doublereal*)
malloc(
sizeof(doublereal)*nn);
761 ipiv = (integer*)
malloc(
sizeof(integer)*
n);
763 for(
i=0;
i<nn;
i++)
a[
i] = _a[
i];
773 cnt += ((
i+1)!=ipiv[
i]);