diff options
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r-- | packages/base/src/Internal/C/lapack-aux.c | 87 |
1 files changed, 87 insertions, 0 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c index 5018a45..dddbefb 100644 --- a/packages/base/src/Internal/C/lapack-aux.c +++ b/packages/base/src/Internal/C/lapack-aux.c | |||
@@ -415,6 +415,93 @@ int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) { | |||
415 | OK | 415 | OK |
416 | } | 416 | } |
417 | 417 | ||
418 | //////////////////// generalized real eigensystem //////////// | ||
419 | |||
420 | int dggev_(char *jobvl, char *jobvr, integer *n, | ||
421 | doublereal *a, integer *lda, doublereal *b, integer *ldb, | ||
422 | doublereal *alphar, doublereal *alphai, doublereal *beta, | ||
423 | doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr, | ||
424 | doublereal *work, | ||
425 | integer *lwork, integer *info); | ||
426 | |||
427 | int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) { | ||
428 | integer n = ar; | ||
429 | REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE); | ||
430 | REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE); | ||
431 | char jobvl = vlp==NULL?'N':'V'; | ||
432 | REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE); | ||
433 | char jobvr = vrp==NULL?'N':'V'; | ||
434 | DEBUGMSG("eig_l_G"); | ||
435 | integer lwork = -1; | ||
436 | integer res; | ||
437 | // ask for optimal lwork | ||
438 | double ans; | ||
439 | dggev_ (&jobvl,&jobvr, | ||
440 | &n, | ||
441 | ap,&n,bp,&n, | ||
442 | (double*)alphap, (double*)alphap+n, betap, | ||
443 | vlp, &n, vrp, &n, | ||
444 | &ans, &lwork, | ||
445 | &res); | ||
446 | lwork = ceil(ans); | ||
447 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
448 | CHECK(!work,MEM); | ||
449 | dggev_ (&jobvl,&jobvr, | ||
450 | &n, | ||
451 | ap,&n,bp,&n, | ||
452 | (double*)alphap, (double*)alphap+n, betap, | ||
453 | vlp, &n, vrp, &n, | ||
454 | work, &lwork, | ||
455 | &res); | ||
456 | CHECK(res,res); | ||
457 | free(work); | ||
458 | OK | ||
459 | } | ||
460 | |||
461 | //////////////////// generalized complex eigensystem //////////// | ||
462 | |||
463 | int zggev_(char *jobvl, char *jobvr, integer *n, | ||
464 | doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, | ||
465 | doublecomplex *alphar, doublecomplex *beta, | ||
466 | doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr, | ||
467 | doublecomplex *work, integer *lwork, | ||
468 | doublereal *rwork, integer *info); | ||
469 | |||
470 | int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) { | ||
471 | integer n = ar; | ||
472 | REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE); | ||
473 | REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE); | ||
474 | char jobvl = vlp==NULL?'N':'V'; | ||
475 | REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE); | ||
476 | char jobvr = vrp==NULL?'N':'V'; | ||
477 | DEBUGMSG("eig_l_GC"); | ||
478 | double *rwork = (double*) malloc(8*n*sizeof(double)); | ||
479 | CHECK(!rwork,MEM); | ||
480 | integer lwork = -1; | ||
481 | integer res; | ||
482 | // ask for optimal lwork | ||
483 | doublecomplex ans; | ||
484 | zggev_ (&jobvl,&jobvr, | ||
485 | &n, | ||
486 | ap,&n,bp,&n, | ||
487 | alphap, betap, | ||
488 | vlp, &n, vrp, &n, | ||
489 | &ans, &lwork, | ||
490 | rwork, &res); | ||
491 | lwork = ceil(ans.r); | ||
492 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
493 | CHECK(!work,MEM); | ||
494 | zggev_ (&jobvl,&jobvr, | ||
495 | &n, | ||
496 | ap,&n,bp,&n, | ||
497 | alphap, betap, | ||
498 | vlp, &n, vrp, &n, | ||
499 | work, &lwork, | ||
500 | rwork, &res); | ||
501 | CHECK(res,res); | ||
502 | free(work); | ||
503 | OK | ||
504 | } | ||
418 | 505 | ||
419 | //////////////////// symmetric real eigensystem //////////// | 506 | //////////////////// symmetric real eigensystem //////////// |
420 | 507 | ||