diff options
Diffstat (limited to 'packages/base/src/Internal/C/lapack-aux.c')
-rw-r--r-- | packages/base/src/Internal/C/lapack-aux.c | 50 |
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, | 357 | int 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 | ||
362 | int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) { | 362 | int 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 * | 405 | int 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 | ||
417 | int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) { | 410 | int 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, | 447 | int 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 | ||
465 | int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) { | 451 | int 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 | 481 | int 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 | ||
503 | int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) { | 485 | int 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, |