summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/C/lapack-aux.c
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Internal/C/lapack-aux.c')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c50
1 files changed, 14 insertions, 36 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 8524962..72c44cb 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -354,12 +354,12 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
354 354
355//////////////////// general complex eigensystem //////////// 355//////////////////// general complex eigensystem ////////////
356 356
357/* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n, 357int zgeev_(char *jobvl, char *jobvr, integer *n,
358 doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl, 358 doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
359 integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work, 359 integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
360 integer *lwork, doublereal *rwork, integer *info); 360 integer *lwork, doublereal *rwork, integer *info);
361 361
362int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) { 362int eig_l_C(OCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
363 integer n = ar; 363 integer n = ar;
364 REQUIRES(ac==n && sn==n, BAD_SIZE); 364 REQUIRES(ac==n && sn==n, BAD_SIZE);
365 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); 365 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -367,18 +367,14 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
367 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE); 367 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
368 char jobvr = vp==NULL?'N':'V'; 368 char jobvr = vp==NULL?'N':'V';
369 DEBUGMSG("eig_l_C"); 369 DEBUGMSG("eig_l_C");
370 doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
371 CHECK(!B,MEM);
372 memcpy(B,ap,n*n*sizeof(doublecomplex));
373 double *rwork = (double*) malloc(2*n*sizeof(double)); 370 double *rwork = (double*) malloc(2*n*sizeof(double));
374 CHECK(!rwork,MEM); 371 CHECK(!rwork,MEM);
375 integer lwork = -1; 372 integer lwork = -1;
376 integer res; 373 integer res;
377 // ask for optimal lwork 374 // ask for optimal lwork
378 doublecomplex ans; 375 doublecomplex ans;
379 //printf("ask zgeev\n");
380 zgeev_ (&jobvl,&jobvr, 376 zgeev_ (&jobvl,&jobvr,
381 &n,B,&n, 377 &n,ap,&n,
382 sp, 378 sp,
383 up,&n, 379 up,&n,
384 vp,&n, 380 vp,&n,
@@ -386,12 +382,10 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
386 rwork, 382 rwork,
387 &res); 383 &res);
388 lwork = ceil(ans.r); 384 lwork = ceil(ans.r);
389 //printf("ans = %d\n",lwork);
390 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); 385 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
391 CHECK(!work,MEM); 386 CHECK(!work,MEM);
392 //printf("zgeev\n");
393 zgeev_ (&jobvl,&jobvr, 387 zgeev_ (&jobvl,&jobvr,
394 &n,B,&n, 388 &n,ap,&n,
395 sp, 389 sp,
396 up,&n, 390 up,&n,
397 vp,&n, 391 vp,&n,
@@ -401,7 +395,6 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
401 CHECK(res,res); 395 CHECK(res,res);
402 free(work); 396 free(work);
403 free(rwork); 397 free(rwork);
404 free(B);
405 OK 398 OK
406} 399}
407 400
@@ -409,12 +402,12 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
409 402
410//////////////////// general real eigensystem //////////// 403//////////////////// general real eigensystem ////////////
411 404
412/* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal * 405int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
413 a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl, 406 a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
414 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work, 407 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
415 integer *lwork, integer *info); 408 integer *lwork, integer *info);
416 409
417int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) { 410int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
418 integer n = ar; 411 integer n = ar;
419 REQUIRES(ac==n && sn==n, BAD_SIZE); 412 REQUIRES(ac==n && sn==n, BAD_SIZE);
420 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); 413 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -422,28 +415,22 @@ int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
422 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE); 415 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
423 char jobvr = vp==NULL?'N':'V'; 416 char jobvr = vp==NULL?'N':'V';
424 DEBUGMSG("eig_l_R"); 417 DEBUGMSG("eig_l_R");
425 double *B = (double*)malloc(n*n*sizeof(double));
426 CHECK(!B,MEM);
427 memcpy(B,ap,n*n*sizeof(double));
428 integer lwork = -1; 418 integer lwork = -1;
429 integer res; 419 integer res;
430 // ask for optimal lwork 420 // ask for optimal lwork
431 double ans; 421 double ans;
432 //printf("ask dgeev\n");
433 dgeev_ (&jobvl,&jobvr, 422 dgeev_ (&jobvl,&jobvr,
434 &n,B,&n, 423 &n,ap,&n,
435 (double*)sp, (double*)sp+n, 424 (double*)sp, (double*)sp+n,
436 up,&n, 425 up,&n,
437 vp,&n, 426 vp,&n,
438 &ans, &lwork, 427 &ans, &lwork,
439 &res); 428 &res);
440 lwork = ceil(ans); 429 lwork = ceil(ans);
441 //printf("ans = %d\n",lwork);
442 double * work = (double*)malloc(lwork*sizeof(double)); 430 double * work = (double*)malloc(lwork*sizeof(double));
443 CHECK(!work,MEM); 431 CHECK(!work,MEM);
444 //printf("dgeev\n");
445 dgeev_ (&jobvl,&jobvr, 432 dgeev_ (&jobvl,&jobvr,
446 &n,B,&n, 433 &n,ap,&n,
447 (double*)sp, (double*)sp+n, 434 (double*)sp, (double*)sp+n,
448 up,&n, 435 up,&n,
449 vp,&n, 436 vp,&n,
@@ -451,37 +438,32 @@ int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
451 &res); 438 &res);
452 CHECK(res,res); 439 CHECK(res,res);
453 free(work); 440 free(work);
454 free(B);
455 OK 441 OK
456} 442}
457 443
458 444
459//////////////////// symmetric real eigensystem //////////// 445//////////////////// symmetric real eigensystem ////////////
460 446
461/* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, 447int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
462 integer *lda, doublereal *w, doublereal *work, integer *lwork, 448 integer *lda, doublereal *w, doublereal *work, integer *lwork,
463 integer *info); 449 integer *info);
464 450
465int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) { 451int eig_l_S(int wantV,DVEC(s),ODMAT(v)) {
466 integer n = ar; 452 integer n = sn;
467 REQUIRES(ac==n && sn==n, BAD_SIZE);
468 REQUIRES(vr==n && vc==n, BAD_SIZE); 453 REQUIRES(vr==n && vc==n, BAD_SIZE);
469 char jobz = wantV?'V':'N'; 454 char jobz = wantV?'V':'N';
470 DEBUGMSG("eig_l_S"); 455 DEBUGMSG("eig_l_S");
471 memcpy(vp,ap,n*n*sizeof(double));
472 integer lwork = -1; 456 integer lwork = -1;
473 char uplo = 'U'; 457 char uplo = 'U';
474 integer res; 458 integer res;
475 // ask for optimal lwork 459 // ask for optimal lwork
476 double ans; 460 double ans;
477 //printf("ask dsyev\n");
478 dsyev_ (&jobz,&uplo, 461 dsyev_ (&jobz,&uplo,
479 &n,vp,&n, 462 &n,vp,&n,
480 sp, 463 sp,
481 &ans, &lwork, 464 &ans, &lwork,
482 &res); 465 &res);
483 lwork = ceil(ans); 466 lwork = ceil(ans);
484 //printf("ans = %d\n",lwork);
485 double * work = (double*)malloc(lwork*sizeof(double)); 467 double * work = (double*)malloc(lwork*sizeof(double));
486 CHECK(!work,MEM); 468 CHECK(!work,MEM);
487 dsyev_ (&jobz,&uplo, 469 dsyev_ (&jobz,&uplo,
@@ -496,17 +478,15 @@ int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) {
496 478
497//////////////////// hermitian complex eigensystem //////////// 479//////////////////// hermitian complex eigensystem ////////////
498 480
499/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex 481int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
500 *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork, 482 *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
501 doublereal *rwork, integer *info); 483 doublereal *rwork, integer *info);
502 484
503int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) { 485int eig_l_H(int wantV,DVEC(s),OCMAT(v)) {
504 integer n = ar; 486 integer n = sn;
505 REQUIRES(ac==n && sn==n, BAD_SIZE);
506 REQUIRES(vr==n && vc==n, BAD_SIZE); 487 REQUIRES(vr==n && vc==n, BAD_SIZE);
507 char jobz = wantV?'V':'N'; 488 char jobz = wantV?'V':'N';
508 DEBUGMSG("eig_l_H"); 489 DEBUGMSG("eig_l_H");
509 memcpy(vp,ap,2*n*n*sizeof(double));
510 double *rwork = (double*) malloc((3*n-2)*sizeof(double)); 490 double *rwork = (double*) malloc((3*n-2)*sizeof(double));
511 CHECK(!rwork,MEM); 491 CHECK(!rwork,MEM);
512 integer lwork = -1; 492 integer lwork = -1;
@@ -514,7 +494,6 @@ int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
514 integer res; 494 integer res;
515 // ask for optimal lwork 495 // ask for optimal lwork
516 doublecomplex ans; 496 doublecomplex ans;
517 //printf("ask zheev\n");
518 zheev_ (&jobz,&uplo, 497 zheev_ (&jobz,&uplo,
519 &n,vp,&n, 498 &n,vp,&n,
520 sp, 499 sp,
@@ -522,7 +501,6 @@ int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
522 rwork, 501 rwork,
523 &res); 502 &res);
524 lwork = ceil(ans.r); 503 lwork = ceil(ans.r);
525 //printf("ans = %d\n",lwork);
526 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); 504 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
527 CHECK(!work,MEM); 505 CHECK(!work,MEM);
528 zheev_ (&jobz,&uplo, 506 zheev_ (&jobz,&uplo,