00001
00002
00003
00004
00005
00006
00007
00008
00016 #include <dims_common.h>
00017 #ifdef __darwin__
00018 typedef int integer;
00019 typedef double doublereal;
00020 #include <stdlib.h>
00021 #else
00022 #include <f2c.h>
00023 #include <malloc.h>
00024 #endif
00025 #include <stdio.h>
00026
00027 #ifdef USE_CLAPACK_INTERFACE
00028 #ifdef __WIN32__
00029 #undef USE_CLAPACK_INTERFACE
00030 #include <blaswrap.h>
00031 #else
00032 #if defined(__darwin__) || defined(__linux__)
00033 extern "C"{
00034 #endif
00035 #include <cblas.h>
00036 #if defined(__darwin__) || defined(__linux__)
00037 }
00038 #endif
00039 #endif
00040 #endif
00041
00042
00043 extern "C"
00044 {
00045 int dgesvx_(
00046 char *fact, char *trans, integer *n,
00047 integer *nrhs, doublereal *a, integer *lda,
00048 doublereal *af, integer *ldaf, integer *ipiv,
00049 char *equed, doublereal *r, doublereal *c,
00050 doublereal *b, integer *ldb, doublereal *x,
00051 integer *ldx, doublereal *rcond, doublereal *ferr,
00052 doublereal *berr, doublereal *work, integer *iwork,
00053 integer *info);
00054
00055 int dgesvd_(
00056 char* jobu, char* jobvt, integer* m, integer* n,
00057 doublereal *a, integer *lda, doublereal *s,
00058 doublereal *u, integer *ldu, doublereal *vt,
00059 integer *ldvt, doublereal* work, integer *lwork,
00060 integer *info);
00061
00062 int dgelss_(
00063 integer* m, integer* n, integer* nrhs,
00064 doublereal* A, integer* lda,
00065 doublereal* B, integer* ldb,
00066 doublereal* S, doublereal* rcond, integer* rank,
00067 doublereal* work, integer* lwork, integer* info);
00068
00069 int dporfs_(
00070 char* uplo, integer* n, integer* nrhs,
00071 doublereal* A, integer* lda,
00072 doublereal* AF, integer* ldaf,
00073 doublereal* b, integer* ldb,
00074 doublereal* x, integer* ldx,
00075 doublereal* ferr, doublereal* berr,
00076 doublereal* work, integer* iwork,
00077 integer* info);
00078
00079 int dpotrf_(char* uplo, integer* n, doublereal* A, integer* lda,
00080 integer* info);
00081
00082 int dpotrs_(char* uplo, integer* n, integer* nrhs,
00083 doublereal* A, integer* lda,
00084 doublereal* x, integer* ldb,
00085 integer* info);
00086
00087 int dposv_(
00088 char* uplo, integer* n, integer* nrhs,
00089 doublereal* A, integer* lda,
00090 doublereal* b, integer* ldb,
00091 integer* info);
00092
00093 int dposvx_(
00094 char* fact, char* uplo, integer* n, integer* nrhs,
00095 doublereal* A, integer* lda, doublereal* AF, integer* ldaf,
00096 char* equed, doublereal* S,
00097 doublereal* B, integer* ldb,
00098 doublereal* X, integer* ldx,
00099 doublereal* rcond, doublereal* ferr, doublereal* berr,
00100 doublereal* work, integer* iwork, integer* info);
00101
00102 int dgesvd_(
00103 char* jobu, char* jobvt, integer* m, integer* n,
00104 doublereal* A, integer* lda,
00105 doublereal* S, doublereal* U, integer* ldu,
00106 doublereal* VT, integer* ldvt,
00107 doublereal* work, integer* lwork, integer* info);
00108
00109 int dsyev_(char* jobz, char* uplo, integer* m, doublereal* a, integer* n, doublereal* w, doublereal* work, integer* lwork, integer* info);
00110
00111 int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *a,
00112 integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
00113 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
00114 integer *lwork, integer *info);
00115
00116 int dgetrf_(integer *m, integer *n, doublereal *a,
00117 integer *lda, integer *ipiv, integer *info);
00118
00119 #ifndef USE_CLAPACK_INTERFACE
00120 int f2c_dscal(
00121 integer* n, doublereal* alpha, doublereal* x, integer* incx);
00122
00123 int f2c_dcopy(
00124 integer* n, doublereal* x, integer* incx, doublereal* y, integer* incy);
00125 double f2c_ddot(
00126 integer* n, doublereal* x, integer* incx, doublereal* y, integer* incy);
00127 int f2c_dgemv(
00128 char* trans, integer* m, integer* n,
00129 doublereal* alpha, doublereal* A, integer* lda,
00130 doublereal* x, integer* incx,
00131 doublereal* beta, doublereal* y, integer* incy);
00132
00133 int f2c_dgemm(
00134 char *transa, char *transb, integer *m, integer *n, integer *k,
00135 doublereal *alpha, doublereal *a, integer *lda,
00136 doublereal *b, integer *ldb,
00137 doublereal *beta, doublereal *c, integer *ldc);
00138
00139 int f2c_dsyrk(
00140 char* uplo, char* trans, integer* n, integer* k,
00141 doublereal* alpha, doublereal* A, integer* lda,
00142 doublereal* beta, doublereal* C, integer* ldc);
00143
00144 int f2c_daxpy(
00145 integer* n, doublereal* alpha, doublereal* x, integer* incx,
00146 doublereal* y, integer* incy);
00147 #endif
00148 }
00149
00150 double dims_dot(double* _x, double* _y, int _n)
00151 {
00152 #ifndef USE_CLAPACK_INTERFACE
00153 integer incx = 1, incy = 1;
00154 integer n = _n;
00155 return f2c_ddot(&n, _x, &incx, _y, &incy);
00156 #else
00157 const int n=_n, incX=1, incY=1;
00158 return cblas_ddot(n, _x, incX, _y, incY);
00159 #endif
00160
00161 }
00162
00163 int dims_copy(double* _x, double* _y, int _n)
00164 {
00165 #ifndef USE_CLAPACK_INTERFACE
00166 integer incx = 1, incy = 1;
00167 integer n = _n;
00168 f2c_dcopy(&n, _x, &incx, _y, &incy);
00169 #else
00170 const int n = _n, incX = 1, incY = 1;
00171 cblas_dcopy(n, _x, incX, _y, incY);
00172 #endif
00173 return 0;
00174 }
00175
00176 int dims_scale_myself(double* _x, double _alpha, int _n)
00177 {
00178 #ifndef USE_CLAPACK_INTERFACE
00179 integer incx = 1;
00180 integer n = _n;
00181 f2c_dscal(&n, &_alpha, _x, &incx);
00182 #else
00183 const int n = _n, incX = 1;
00184 const double alpha = _alpha;
00185 cblas_dscal(n, alpha, _x, incX);
00186 #endif
00187 return 0;
00188 }
00189
00190 int dims_scale(double* _x, double _alpha, int _n, double*_y)
00191 {
00192 #ifndef USE_CLAPACK_INTERFACE
00193 integer incx = 1, incy = 1;
00194 integer n = _n;
00195 f2c_dcopy(&n, _x, &incx, _y, &incy);
00196 f2c_dscal(&n, &_alpha, _y, &incx);
00197 #else
00198 const int incX = 1, incY = 1;
00199 cblas_dcopy(_n, _x, incX, _y, incY);
00200 cblas_dscal(_n, _alpha, _x, incX);
00201 #endif
00202 return 0;
00203 }
00204
00205 int dims_dgemv(double* _A, int _m, int _n, double* _x, double* _y)
00206 {
00207 #ifndef USE_CLAPACK_INTERFACE
00208 char trans = 'N';
00209 integer m = _m;
00210 integer n = _n;
00211 integer lda = _m;
00212 doublereal alpha = 1.0;
00213 integer incx = 1;
00214 doublereal beta = 0.0;
00215 integer incy = 1;
00216 f2c_dgemv(&trans, &m, &n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
00217 #else
00218 cblas_dgemv(CblasColMajor, CblasNoTrans,
00219 _m, _n, 1.0, _A,_m, _x, 1, 0.0, _y, 1);
00220 #endif
00221 return 0;
00222 }
00223
00224 int dims_dgemv_tran(double* _A, int _m, int _n, double* _x, double* _y)
00225 {
00226 #ifndef USE_CLAPACK_INTERFACE
00227 char trans = 'T';
00228 integer m = _m;
00229 integer n = _n;
00230 integer lda = _m;
00231 doublereal alpha = 1.0;
00232 integer incx = 1;
00233 doublereal beta = 0.0;
00234 integer incy = 1;
00235 f2c_dgemv(&trans, &m, &n, &alpha, _A, &lda, _x, &incx, &beta, _y, &incy);
00236 #else
00237 cblas_dgemv(CblasColMajor, CblasTrans,
00238 _m, _n, 1.0, _A, _m, _x, 1, 0.0, _y, 1);
00239 #endif
00240 return 0;
00241 }
00242
00243 int dims_dgemm(double* _A, double* _B, int _m, int _n, int _k, double* _C)
00244 {
00245 #ifndef USE_CLAPACK_INTERFACE
00246 char transa = 'N';
00247 char transb = 'N';
00248 integer m = _m;
00249 integer n = _n;
00250 integer k = _k;
00251 doublereal alpha = 1.0;
00252 integer lda = _m;
00253 integer ldb = _k;
00254 doublereal beta = 0.0;
00255 integer ldc = _m;
00256 f2c_dgemm(&transa, &transb, &m, &n, &k, &alpha, _A, &lda, _B, &ldb, &beta, _C, &ldc);
00257 #else
00258 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
00259 _m, _n, _k, 1.0,
00260 _A, _m,
00261 _B, _k,
00262 0.0,
00263 _C, _m
00264 );
00265 #endif
00266 return 0;
00267 }
00268
00269 int dims_dsyrk(double* _A, int _n, int _k, double* _C)
00270 {
00271 #ifndef USE_CLAPACK_INTERFACE
00272 char uplo = 'U';
00273 char trans = 'N';
00274 integer n = _n;
00275 integer k = _k;
00276 doublereal alpha = 1.0;
00277 integer lda = _n;
00278 doublereal beta = 0.0;
00279 integer ldc = _n;
00280
00281 f2c_dsyrk(&uplo, &trans, &n, &k, &alpha, _A, &lda, &beta, _C, &ldc);
00282 #else
00283 cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans,
00284 _n, _k, 1.0, _A, _n, 0.0, _C, _n);
00285 #endif
00286 return 0;
00287 }
00288
00289 int dims_dsyrk_trans_first(double* _A, int _n, int _k, double* _C)
00290 {
00291 #ifndef USE_CLAPACK_INTERFACE
00292 char uplo = 'U';
00293 char trans = 'T';
00294 integer n = _n;
00295 integer k = _k;
00296 doublereal alpha = 1.0;
00297 integer lda = _k;
00298 doublereal beta = 0.0;
00299 integer ldc = _n;
00300 f2c_dsyrk(&uplo, &trans, &n, &k, &alpha, _A, &lda, &beta, _C, &ldc);
00301 #else
00302 cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans,
00303 _n, _k, 1.0, _A, _n, 0.0, _C, _n);
00304 #endif
00305 return 0;
00306 }
00307
00308 int dims_daxpy(int _n, double _alpha, double* _x, double* _y)
00309 {
00310 #ifndef USE_CLAPACK_INTERFACE
00311 integer n = _n;
00312 doublereal alpha = _alpha;
00313 integer incx = 1;
00314 integer incy = 1;
00315 f2c_daxpy(&n, &alpha, _x, &incx, _y, &incy);
00316 #else
00317 cblas_daxpy(_n, _alpha, _x, 1, _y, 1);
00318 #endif
00319 return 0;
00320 }
00321
00322 int dims_dporfs(double* _a, double* _x, double* _b, int _m, int _nrhs)
00323 {
00324 char uplo = 'U';
00325 integer n = _m;
00326 integer nrhs = _nrhs;
00327 integer lda = _m;
00328 integer ldaf = _m;
00329 integer ldb = _m;
00330 integer ldx = _m;
00331 integer info;
00332 integer a_ndata = n*n;
00333 integer b_ndata = n*nrhs;
00334 integer i;
00335 doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
00336 doublereal* u = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
00337 doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
00338 doublereal* x = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
00339 doublereal* ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00340 doublereal* berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00341 doublereal* work = (doublereal *)malloc(sizeof(doublereal)*3*n);
00342 integer* iwork = (integer*)malloc(sizeof(integer)*n);
00343
00344 for(i=0; i<a_ndata; i++)
00345 {
00346 a[i] = _a[i];
00347 u[i] = _a[i];
00348 }
00349 for(i=0; i<b_ndata; i++)
00350 {
00351 b[i] = _b[i];
00352 x[i] = _b[i];
00353 }
00354 dpotrf_(&uplo, &n, u, &lda, &info);
00355 dpotrs_(&uplo, &n, &nrhs, u, &lda, x, &ldb, &info);
00356 #if 0
00357 printf("info=%d\n", info);
00358 printf("old x = ");
00359 for(i=0; i<b_ndata; i++)
00360 {
00361 printf("\t%.8g", x[i]);
00362 }
00363 printf("\n");
00364 #endif
00365
00366 dporfs_(&uplo, &n, &nrhs, a, &lda, u, &ldaf, b, &ldb, x, &ldx, ferr, berr, work, iwork, &info);
00367 for(i=0; i<b_ndata; i++)
00368 {
00369 _x[i] = x[i];
00370 }
00371 #if 0
00372 printf("info=%d\n", info);
00373 printf("new x = ");
00374 for(i=0; i<b_ndata; i++)
00375 {
00376 printf("\t%.8g", x[i]);
00377 }
00378 printf("\n");
00379 #endif
00380 free(a);
00381 free(u);
00382 free(b);
00383 free(x);
00384 free(ferr);
00385 free(berr);
00386 free(work);
00387 free(iwork);
00388 return info;
00389 }
00390
00391 int dims_dposvx(double* _a, double* _x, double* _b, int _m, int _nrhs, double* _rcond)
00392 {
00393 char fact = 'E';
00394 char uplo = 'U';
00395 integer n = _m;
00396 integer nrhs = _nrhs;
00397 integer lda = _m;
00398 integer ldaf = _m;
00399 char equed = 'N';
00400 integer ldb = _m;
00401 integer ldx = _m;
00402 doublereal rcond;
00403 integer info;
00404 integer a_ndata = n*n;
00405 integer b_ndata = n*nrhs;
00406 integer i;
00407 doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
00408 doublereal* af = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
00409 doublereal* s = (doublereal *)malloc(sizeof(doublereal)*n);
00410 doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
00411 doublereal* x = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
00412 doublereal* ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00413 doublereal* berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00414 doublereal* work = (doublereal *)malloc(sizeof(doublereal)*3*n);
00415 integer* iwork = (integer*)malloc(sizeof(integer)*n);
00416
00417 for(i=0; i<a_ndata; i++)
00418 {
00419 a[i] = _a[i];
00420 }
00421 for(i=0; i<b_ndata; i++)
00422 {
00423 b[i] = _b[i];
00424 }
00425 dposvx_(&fact, &uplo, &n, &nrhs, a, &lda, af, &ldaf, &equed, s, b, &ldb, x, &ldx, &rcond, ferr, berr, work, iwork, &info);
00426 for(i=0; i<b_ndata; i++)
00427 {
00428 _x[i] = x[i];
00429 }
00430 #if 0
00431 printf("scale=[");
00432 for(i=0; i<n; i++)
00433 {
00434 printf("\t%.8g", s[i]);
00435 }
00436 printf("]\n");
00437 #endif
00438 *_rcond = rcond;
00439 free(a);
00440 free(af);
00441 free(s);
00442 free(b);
00443 free(x);
00444 free(ferr);
00445 free(berr);
00446 free(work);
00447 free(iwork);
00448 return info;
00449 }
00450
00451 int dims_dposv(double* _a, double* _x, double* _b, int _m, int _nrhs)
00452 {
00453 char uplo = 'U';
00454 integer n = _m;
00455 integer nrhs = _nrhs;
00456 integer lda = _m;
00457 integer ldb = _m;
00458 integer info;
00459 integer a_ndata = n*n;
00460 integer b_ndata = n*nrhs;
00461 integer i;
00462 doublereal* a = (doublereal *)malloc(sizeof(doublereal)*a_ndata);
00463 doublereal* b = (doublereal *)malloc(sizeof(doublereal)*b_ndata);
00464
00465 for(i=0; i<a_ndata; i++)
00466 {
00467 a[i] = _a[i];
00468 }
00469 for(i=0; i<b_ndata; i++)
00470 {
00471 b[i] = _b[i];
00472 }
00473 dposv_(&uplo, &n, &nrhs, a, &lda, b, &ldb, &info);
00474 for(i=0; i<b_ndata; i++)
00475 {
00476 _x[i] = b[i];
00477 }
00478 free(a);
00479 free(b);
00480 return info;
00481 }
00482
00483 int dims_dgesvx(double* _a, double* _x, double* _b, int _n, int _nrhs)
00484 {
00485 char fact = 'N';
00486 char trans = 'N';
00487 integer n = (integer)_n;
00488 integer nrhs = (integer)_nrhs;
00489 doublereal* a = (doublereal*)_a;
00490 integer lda = n;
00491 doublereal* af = (doublereal*)malloc(sizeof(doublereal)*n*n);
00492 integer ldaf = n;
00493 integer *ipiv=(integer *)malloc(sizeof(integer)*n);
00494 char equed = 'N';
00495 doublereal *r=(doublereal *)malloc(sizeof(doublereal)*n);
00496 doublereal *c=(doublereal *)malloc(sizeof(doublereal)*n);
00497 doublereal* b = (doublereal*)_b;
00498 integer ldb = n;
00499 doublereal *x = (doublereal *)_x;
00500 integer ldx = n;
00501
00502 doublereal rcond;
00503 doublereal *ferr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00504 doublereal *berr = (doublereal *)malloc(sizeof(doublereal)*nrhs);
00505 doublereal *work = (doublereal *)malloc(sizeof(doublereal)*4*n);
00506 integer *iwork = (integer *)malloc(sizeof(integer)*n);
00507
00508 integer info;
00509 dgesvx_(&fact, &trans, &n, &nrhs, a, &lda, af, &ldaf, ipiv,
00510 &equed, r, c, b, &ldb, x, &ldx, &rcond, ferr, berr,
00511 work, iwork, &info);
00512 free(ipiv);
00513 free(r);
00514 free(c);
00515 free(iwork);
00516 free(work);
00517 free(berr);
00518 free(ferr);
00519 free(af);
00520 return info;
00521 }
00522
00523
00524 int dims_dgelss(double* _a, double* _x, double* _b, int _m, int _n, int _nrhs,
00525 double* _s, int* _rank, int _lwork)
00526 {
00527 int i, mm, tmp1, tmp2;
00528 integer m = (integer)_m;
00529 integer n = (integer)_n;
00530 integer nrhs = (integer)_nrhs;
00531 integer lda, ldb, rank, lwork, info;
00532 doublereal *a, *b, *s, rcond, *work;
00533 mm = _m*_n;
00534 a = (doublereal*)malloc(sizeof(doublereal)*mm);
00535 for(i=0; i<mm; i++) a[i] = _a[i];
00536 lda = MAX(1, _m);
00537 ldb = MAX3(1, _m, _n);
00538 mm = ldb*_nrhs;
00539 b = (doublereal*)malloc(sizeof(doublereal)*mm);
00540 for(i=0; i<mm; i++) b[i] = _b[i];
00541 s = (doublereal*)_s;
00542 rcond = -1;
00543 if(_lwork > 0) lwork = _lwork;
00544 else
00545 {
00546 lwork = MIN(_m, _n);
00547 lwork *= 3;
00548 tmp1 = MIN(_m, _n);
00549 tmp2 = MAX(_m, _n);
00550 tmp1 *= 2;
00551 lwork += MAX3(tmp1, tmp2, _nrhs);
00552 lwork *= 5;
00553 }
00554 work = (doublereal*)malloc(sizeof(doublereal)*lwork);
00555 dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank,
00556 work, &lwork, &info);
00557
00558 mm = _n*_nrhs;
00559 for(i=0; i<mm; i++) _x[i] = b[i];
00560 *_rank = (int)rank;
00561
00562 free(a);
00563 free(b);
00564 free(work);
00565 return info;
00566 }
00567
00568 int dims_svd(double* _a, int _m, int _n, double* _u, double* _sigma, double* _vt)
00569 {
00570 char jobu = 'A';
00571 char jobvt = 'A';
00572 int mn = _m * _n, i, large = MAX(_m, _n), small = MIN(_m, _n);
00573 integer m = (integer)_m;
00574 integer n = (integer)_n;
00575 integer lda = m;
00576 integer ldu = m;
00577 integer ldvt = n;
00578 integer lwork, info;
00579 doublereal* a, *work;
00580 doublereal* u = (doublereal*)_u;
00581 doublereal* s = (doublereal*)_sigma;
00582 doublereal* vt = (doublereal*)_vt;
00583
00584 a = (doublereal*)malloc(sizeof(doublereal) * mn);
00585 for(i=0; i<mn; i++) a[i] = _a[i];
00586
00587 lwork = MAX(3*small+large, 5*small-4);
00588 lwork *= 5;
00589
00590 work = (doublereal*)malloc(sizeof(doublereal) * lwork);
00591
00592 dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);
00593
00594 free(a);
00595 free(work);
00596 return info;
00597 }
00598
00599
00600 int dims_eigs(int _n, double *_a, double *w)
00601 {
00602 char jobz='V';
00603 char uplo='U';
00604 integer lwork;
00605 double *work;
00606 double *a;
00607 integer info;
00608 int i;
00609 integer n = (integer)_n;
00610
00611 a=(double*)malloc(sizeof(double)*n*n);
00612 lwork=3*n;
00613 lwork=lwork*5;
00614 work=(double*)malloc(sizeof(double)*lwork);
00615
00616 for(i=0;i<n*n;i++){
00617 a[i]=_a[i];
00618 }
00619 dsyev_(&jobz,&uplo,&n,a,&n,w,work,&lwork,&info);
00620 for(i=0;i<n*n;i++)
00621 {
00622 _a[i]=a[i];
00623 }
00624 free(a);
00625 free(work);
00626 return info;
00627 }
00628
00629 int dims_eigs2(int _n, double *_a, double *w)
00630 {
00631 char jobz='V';
00632 char uplo='U';
00633 integer lwork;
00634 double *work=NULL;
00635 double *a=NULL;
00636 integer info;
00637 int i;
00638 integer n = (integer)_n;
00639
00640 if(!a)
00641 a=(double*)malloc(sizeof(double)*n*n);
00642 lwork=3*n;
00643 lwork=lwork*5;
00644 if(!work)
00645 work=(double*)malloc(sizeof(double)*lwork);
00646
00647 for(i=0;i<n*n;i++)
00648 {
00649 a[i]=_a[i];
00650 }
00651 dsyev_(&jobz,&uplo,&n,a,&n,w,work,&lwork,&info);
00652 for(i=0;i<n*n;i++)
00653 {
00654 _a[i]=a[i];
00655 }
00656 free(a);
00657 free(work);
00658
00659 return info;
00660 }
00661
00662 int dims_dgeev(int _n, double* _a, double* _wr, double* _wi, double* _vr)
00663 {
00664 int i, nn;
00665 char jobvl = 'N';
00666 char jobvr = 'V';
00667 integer n, lda, ldvl, ldvr, lwork, info;
00668 doublereal *a, *wr, *wi, *vl, *vr, *work;
00669
00670 n = (integer)_n;
00671 lda = n;
00672 ldvl = 1;
00673 ldvr = n;
00674 lwork = 4*n;
00675 nn = _n*_n;
00676
00677 a = (doublereal*)malloc(sizeof(doublereal)*nn);
00678 wr = (doublereal*)malloc(sizeof(doublereal)*n);
00679 wi = (doublereal*)malloc(sizeof(doublereal)*n);
00680 vl = (doublereal*)malloc(sizeof(doublereal)*ldvl*n);
00681 vr = (doublereal*)malloc(sizeof(doublereal)*ldvr*n);
00682 work = (doublereal *)malloc(sizeof(doublereal)*lwork);
00683
00684 for(i=0; i<nn; i++) a[i] = _a[i];
00685
00686 dgeev_(&jobvl, &jobvr, &n,
00687 a, &lda,
00688 wr, wi,
00689 vl, &ldvl,
00690 vr, &ldvr,
00691 work, &lwork, &info);
00692
00693 for(i=0; i<n; i++) _wr[i] = wr[i];
00694 for(i=0; i<n; i++) _wi[i] = wi[i];
00695 for(i=0; i<nn; i++) _vr[i] = vr[i];
00696 free(a);
00697 free(wr);
00698 free(wi);
00699 free(vl);
00700 free(vr);
00701 free(work);
00702 return info;
00703 }
00704
00705 int dims_dgeev_simple(int _n, double* _a, double* _wr, double* _wi)
00706 {
00707 int i, nn;
00708 char jobvl = 'N';
00709 char jobvr = 'N';
00710 integer n, lda, ldvl, ldvr, lwork, info;
00711 doublereal *a, *wr, *wi, *vl, *vr, *work;
00712
00713 n = (integer)_n;
00714 lda = n;
00715 ldvl = 1;
00716 ldvr = n;
00717 lwork = 3*n;
00718 nn = _n*_n;
00719
00720 a = (doublereal*)malloc(sizeof(doublereal)*nn);
00721 wr = (doublereal*)malloc(sizeof(doublereal)*n);
00722 wi = (doublereal*)malloc(sizeof(doublereal)*n);
00723 vl = (doublereal*)malloc(sizeof(doublereal)*ldvl*n);
00724 vr = (doublereal*)malloc(sizeof(doublereal)*ldvr*n);
00725 work = (doublereal *)malloc(sizeof(doublereal)*lwork);
00726
00727 for(i=0; i<nn; i++) a[i] = _a[i];
00728
00729 dgeev_(&jobvl, &jobvr, &n,
00730 a, &lda,
00731 wr, wi,
00732 vl, &ldvl,
00733 vr, &ldvr,
00734 work, &lwork, &info);
00735
00736 for(i=0; i<n; i++) _wr[i] = wr[i];
00737 for(i=0; i<n; i++) _wi[i] = wi[i];
00738
00739 free(a);
00740 free(wr);
00741 free(wi);
00742 free(vl);
00743 free(vr);
00744 free(work);
00745 return info;
00746 }
00747
00748 int dims_det(int _n, double* _a, double* _x)
00749 {
00750 int i, nn, cnt = 0;
00751 const double sign[2] = {1.0, -1.0};
00752
00753 integer n, info;
00754 doublereal *a;
00755 integer* ipiv;
00756
00757 n = (integer)_n;
00758 nn = _n*_n;
00759
00760 a = (doublereal*)malloc(sizeof(doublereal)*nn);
00761 ipiv = (integer*)malloc(sizeof(integer)*n);
00762
00763 for(i=0; i<nn; i++) a[i] = _a[i];
00764
00765 dgetrf_(&n, &n,
00766 a, &n,
00767 ipiv, &info);
00768
00769 *_x = 1.0;
00770 for(i=0; i<n; i++)
00771 {
00772 *_x *= a[i*n+i];
00773 cnt += ((i+1)!=ipiv[i]);
00774 }
00775 *_x *= sign[cnt%2];
00776
00777 free(a);
00778 free(ipiv);
00779 return info;
00780 }