, 11 min read
Computing Generalized Eigenvalues using GUPTRI
This post shows how to call the GUPTRI algorithm from C. The GUPTRI algorithm solves the generalized eigenvalue problem $A x = \lambda B x$. The algorithm can properly handle the cases when $A$, $B$, or both matrices are singular.
Bibliography.
- Demmel, James and Kagstrom, Bo: “The Generalized Schur Decomposition of an Arbitrary Pencil $A – \lambda B$: Robust Software with Error Bounds and Applications. Part I: theory and algorithms”, PDF, ACM TOMS, Vol. 19, Issue 2, June 1993, pp.160–174
- Demmel, James and Kagstrom, Bo: “The Generalized Schur Decomposition of an Arbitrary Pencil $A – \lambda B$: Robust Software with Error Bounds and Applications. Part II: Software and Applications”, PDF, ACM TOMS, Vol. 19, No. 2, June 1993, pp.175–201
1. QZ algorithm fails for singular matrix pencils
octave:4> m0=[-653.0 / 200, 4203.0 / 25400; -1.63, 10527.0 / 127000]
m0 =
-3.265000 0.165472
-1.630000 0.082890
octave:5> m1=[109.0 / 200, 94929.0 / 127000; 0, 0.387]
m1 =
0.5450 0.7475
0 0.3870
octave:6> eig(m0,-m1)
ans =
6.586532109869343e-02
-6.586532109869245e-02
But embedding this 2x2 problem into a 3x3 problem gives entirely wrong results:
octave:9> m03=[0,0,0; 0,-653.0 / 200, 4203.0 / 25400; 0,-1.63, 10527.0 / 127000]
m03 =
0 0 0
0 -3.265000000000000 0.165472440944882
0 -1.630000000000000 0.082889763779528
octave:10> m13=[0,0,0; 0,109.0 / 200, 94929.0 / 127000; 0,0, 0.387]
m13 =
0 0 0
0 0.545000000000000 0.747472440944882
0 0 0.387000000000000
octave:11> eig(m03,-m13)
ans =
NaN
5.990825688073394
NaN
Also, if you ask Grok or Kimi about these kind of problems, these AI chatbots employ the QZ algorithm and therefore give wrong answers.
2. Compiling GUPTRI from source
Download GUPTRI from Netlib here: https://www.netlib.org/linalg/guptri
The file is a shell archive.
After unpacking some of the Fortran routines need to be compiled with the flag -std=legacy.
gfortran -c -O2 zblas.f
gfortran -c -O2 zbnd.f
gfortran -c -O2 zcmatmlr.f
gfortran -c -O2 zftest1.f
gfortran -std=legacy -c -O2 zgschurm.f
gfortran -c -O2 zguptri.f
gfortran -c -O2 zlinpack.f
gfortran -c -O2 zlistr.f
gfortran -std=legacy -c -O2 zmiscl.f
gfortran -std=legacy -c -O2 zqz.f
gfortran -c -O2 zrcsvdc.f
gfortran -c -O2 zreorder.f
gfortran -c -O2 zrzstr.f
Above routines must be linked to your calling program.
3. Testing GUPTRI from C
The actual GUPTRI source code is in Fortran.
Compiling goes like this:
cc -Wall guptritst.c $G/zguptri.o $G/zrzstr.o $G/zcmatmlr.o $G/zqz.o $G/zmiscl.o $G/zlinpack.o $G/zblas.o $G/zrcsvdc.o $G/zlistr.o -o guptritst -lgfortran -lm
$G is the location of the GUPTRI object files.
We call the Fortran GUPTRI algorithm this from C.
/* Test program for using GUPTR from Demmel + Kagstrom
The Generalized Schur Decomposition of an Arbitrary Pencil A - lambda B:
Robust Software with Error Bounds and Applications.
Part II: Software and Applications
ACM TOMS, Vol. 19, No. 2, June 1993, pp.175--201
*/
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <complex.h>
extern void guptri_ (complex a[], complex b[], int *ldab, int *m, int *n, double *epsu, double *gap, int *zero,
complex p[], int *ldp, complex q[], int *ldq,
double *adelta, double *bdelta, int *rtre, int *rtce, int *zrre, int *zrce,
int *fnre, int *fnce, int *inre, int *ince, int pstruc[], int struc[],
complex work[], int kstr[], int *info);
Below struct is required if you want to access the COMMON block variables in Fortran.
extern struct {
int idbg[20];
int outunit;
} debug2_;
We test with a classical eigenvalue problem from a Numerical Analysis book.
#define N 6
#define NSQ 36
/* Generalized eigenvalues are
( 12.316876 0.000000)
( 9.016136 0.000000)
( 6.000000 0.000000)
( 2.983864 0.000000)
( -0.316876 0.000000)
*/
double stoer[] = { // Stoer/Burlisch (1970):``Introduction to Numerical Analysis", page 374
12, 1, 0, 0, 0,
1, 9, 1, 0, 0,
0, 1, 6, 1, 0,
0, 0, 1, 3, 1,
0, 0, 0, 1, 0
};
And, of course, we want to solve the Mihelcic4 problem.
double mihelcic4A2[] = {
-653.0 / 200, 4203.0 / 25400,
-1.63, 10527.0 / 127000
};
double mihelcic4B2[] = {
109.0 / 200, 94929.0 / 127000,
0, 0.387
};
Some helper function.
void zmatprt (const char s[], int n, complex a[]) {
int i, j;
puts(s);
for (i=0; i<n; ++i) {
for (j=0; j<n; ++j)
printf(" (%10.4f, %10.4f)",creal(a[i+j*n]),cimag(a[i+j*n]));
puts("");
}
}
Below function uses the logic detailed in zguptri.f around label 2005.
void zeigvprt (int n, complex a[], complex b[], int rtre, int rtce, int zrre, int zrce, int fnre, int fnce, int inre) {
int i, j;
complex z;
int djordz = zrre - rtre;
int djordi = inre - fnre;
int dimreg = fnre - zrre;
int allreg = dimreg + djordz + djordi;
int rowb = zrre - djordz;
int colb = zrce - djordz;
int rowe = rowb + allreg;
printf("eigenvalues before reordering:\n");
for (i=rowb; i<rowe; ++i) {
j = colb + i - rowb;
if (cabs(b[i+j*n]) == 0.0) {
printf("\tinfinite eigenvalue: (%g,%g), (%g,%g)\n",creal(a[i+j*n]),cimag(a[i+j*n]), creal(b[i+j*n]),cimag(b[i+j*n]));
} else {
z = a[i+j*n] / b[i+j*n];
printf("\teigenvalue = (%g, %g)\n",creal(z),cimag(z));
}
}
}
Finally, the main routine.
int main (int argc, char *argv[]) {
complex a[NSQ], b[NSQ], p[NSQ], q[NSQ], work[2*N+NSQ+NSQ+6*N+N+10];
double epsu=100*DBL_EPSILON, gap=1000, adelta, bdelta;
int ldab=N, m=N, n=N, zero=0, ldp=N, ldq=N, rtre, rtce, zrre, zrce, fnre, fnce, inre, ince, info;
int pstruc[4], struc[N], kstr[4*N];
int i, j;
ldab = 5, m = 5, n = 5;
for (i=0; i<n; ++i) {
for (j=0; j<n; ++j)
a[i+j*n] = stoer[i*n+j], // transpose
b[i+j*n] = 0;
b[i+i*n] = 1;
}
debug2_.idbg[1] = 1;
debug2_.outunit = 6;
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zmatprt("stoer output",N,a);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
ldab = 2, m = 2, n = 2;
for (i=0; i<n; ++i)
for (j=0; j<n; ++j)
a[i+j*n] = mihelcic4A2[i*n+j], // transpose
b[i+j*n] = 0 - mihelcic4B2[i*n+j]; // dto.
zmatprt("\n\n* * *\nMihelcic4 2x2 input A",n,a);
zmatprt("Mihelcic4 2x2 input B",n,b);
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
zmatprt("Mihelcic4 2x2 output",n,a);
zmatprt("Mihelcic4 2x2 output",n,b);
ldab = 3, m = 3, n = 3;
memset(a,0,n*n*sizeof(complex));
memset(b,0,n*n*sizeof(complex));
for (i=1; i<n; ++i)
for (j=1; j<n; ++j)
a[i+j*n] = mihelcic4A2[(i-1)*2+(j-1)], // transpose
b[i+j*n] = 0 - mihelcic4B2[(i-1)*2+(j-1)]; // dto.
zmatprt("\n\n* * *\nMihelcic4 3x3 input A",n,a);
zmatprt("Mihelcic4 3x3 input B",n,b);
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
zmatprt("Mihelcic4 3x3 output",n,a);
zmatprt("Mihelcic4 3x3 output",n,b);
ldab = 4, m = 4, n = 4;
memset(a,0,n*n*sizeof(complex));
memset(b,0,n*n*sizeof(complex));
for (i=1; i<3; ++i)
for (j=1; j<3; ++j)
a[i+j*n] = mihelcic4A2[(i-1)*2+(j-1)], // transpose
b[i+j*n] = 0 - mihelcic4B2[(i-1)*2+(j-1)]; // dto.
zmatprt("\n\n* * *\nMihelcic4 4x4 input A",n,a);
zmatprt("Mihelcic4 4x4 input B",n,b);
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
zmatprt("Mihelcic4 4x4 output",n,a);
zmatprt("Mihelcic4 4x4 output",n,b);
ldab = 5, m = 5, n = 5;
memset(a,0,n*n*sizeof(complex));
memset(b,0,n*n*sizeof(complex));
for (i=1; i<3; ++i)
for (j=1; j<3; ++j)
a[i+j*n] = mihelcic4A2[(i-1)*2+(j-1)], // transpose
b[i+j*n] = 0 - mihelcic4B2[(i-1)*2+(j-1)]; // dto.
zmatprt("\n\n* * *\nMihelcic4 5x5 input A",n,a);
zmatprt("Mihelcic4 5x5 input B",n,b);
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
zmatprt("Mihelcic4 5x5 output",n,a);
zmatprt("Mihelcic4 5x5 output",n,b);
ldab = 2, m = 2, n = 2;
a[0]=1, a[2]=0,
a[1]=0, a[3]=0;
b[0]=1, b[2]=0,
b[1]=0, b[3]=1;
guptri_(a,b,&ldab,&m,&n,&epsu,&gap,&zero,p,&ldp,q,&ldq,&adelta,&bdelta,
&rtre,&rtce,&zrre,&zrce,&fnre,&fnce,&inre,&ince,pstruc,struc,work,kstr,&info);
printf("guptri_: info=%d\n",info);
printf("guptri_: rtre=%d, rtce=%d, zrre=%d, zrce=%d, fnre=%d, fnce=%d, inre=%d, ince=%d\n",
rtre,rtce,zrre,zrce,fnre,fnce,inre,ince);
printf("guptri_: dim regular part = %d, zero mult = %d, infty mult = %d, total mult = %d\n",
inre-rtre,zrre-rtre,inre-fnre,fnre-zrre);
zeigvprt(n,a,b,rtre,rtce,zrre,zrce,fnre,fnce,inre);
return 0;
}
The output is like this.
Stoer matrix.
guptri_: info=0
guptri_: rtre=0, rtce=0, zrre=0, zrce=0, fnre=5, fnce=5, inre=5, ince=5
guptri_: dim regular part = 5, zero mult = 0, infty mult = 0, total mult = 5
stoer output
( 12.3169, 0.0000) ( 9.0161, 0.0000) ( 6.0000, 0.0000) ( 2.9839, 0.0000) ( -0.3169, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( -0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
eigenvalues before reordering:
eigenvalue = (12.3169, 0)
eigenvalue = (9.01614, 0)
eigenvalue = (6, 0)
eigenvalue = (2.98386, 0)
eigenvalue = (-0.316876, 0)
Mihelcic4 as 2x2 matrix.
* * *
Mihelcic4 2x2 input A
( -3.2650, 0.0000) ( 0.1655, 0.0000)
( -1.6300, 0.0000) ( 0.0829, 0.0000)
Mihelcic4 2x2 input B
( -0.5450, 0.0000) ( -0.7475, 0.0000)
( 0.0000, 0.0000) ( -0.3870, 0.0000)
guptri_: info=0
guptri_: rtre=0, rtce=0, zrre=0, zrce=0, fnre=2, fnce=2, inre=2, ince=2
guptri_: dim regular part = 2, zero mult = 0, infty mult = 0, total mult = 2
eigenvalues before reordering:
eigenvalue = (0.0658653, 0)
eigenvalue = (-0.0658653, 0)
Mihelcic4 2x2 output
( 0.0574, 0.0000) ( -3.6535, -0.0000)
( 0.0000, -0.0000) ( -0.0159, 0.0000)
Mihelcic4 2x2 output
( 0.8721, 0.0000) ( -0.4318, -0.0000)
( 0.0000, 0.0000) ( 0.2418, 0.0000)
Mihelcic4 as 3x3 matrix.
* * *
Mihelcic4 3x3 input A
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -3.2650, 0.0000) ( 0.1655, 0.0000)
( 0.0000, 0.0000) ( -1.6300, 0.0000) ( 0.0829, 0.0000)
Mihelcic4 3x3 input B
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -0.5450, 0.0000) ( -0.7475, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.3870, 0.0000)
guptri_: info=0
guptri_: rtre=0, rtce=1, zrre=0, zrce=1, fnre=2, fnce=3, inre=2, ince=3
guptri_: dim regular part = 2, zero mult = 0, infty mult = 0, total mult = 2
eigenvalues before reordering:
eigenvalue = (-0.0658653, 0)
eigenvalue = (0.0658653, 0)
Mihelcic4 3x3 output
( 0.0000, 0.0000) ( -0.0565, -0.0000) ( -3.6535, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0162, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 3x3 output
( 0.0000, 0.0000) ( 0.8583, 0.0000) ( -0.4566, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.2457, 0.0000)
( 0.0000, 0.0000) ( -0.0000, 0.0000) ( -0.0000, 0.0000)
Mihelcic4 as 4x4 matrix.
* * *
Mihelcic4 4x4 input A
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -3.2650, 0.0000) ( 0.1655, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -1.6300, 0.0000) ( 0.0829, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 4x4 input B
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -0.5450, 0.0000) ( -0.7475, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.3870, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
guptri_: info=0
guptri_: rtre=0, rtce=2, zrre=0, zrce=2, fnre=2, fnce=4, inre=2, ince=4
guptri_: dim regular part = 2, zero mult = 0, infty mult = 0, total mult = 2
eigenvalues before reordering:
eigenvalue = (-0.0658653, 0)
eigenvalue = (0.0658653, 0)
Mihelcic4 4x4 output
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.0565, -0.0000) ( -3.6535, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0162, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 4x4 output
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.8583, 0.0000) ( -0.4566, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.2457, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.0000, 0.0000) ( -0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 as 5x5 matrix.
* * *
Mihelcic4 5x5 input A
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -3.2650, 0.0000) ( 0.1655, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -1.6300, 0.0000) ( 0.0829, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 5x5 input B
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( -0.5450, 0.0000) ( -0.7475, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.3870, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
guptri - workspace for rzstr - 5 5 5 5 1
6 31 37 67 42 47 92 52 57 62
guptri - m,n,epsu= 5 5 0.222045D-13
reduction 1
kstr, last= 1
1
0
0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 2
1 2
-1 0
-1 0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
final kstr=
1 2 3 4
-1 -1 5 -1
-1 -1 5 -1
nsumrz= 0
rsumrz= 0
djordz= 0
nsumli= 0
rsumli= 0
djordi= 0
dimreg= 5
final pstruc= 0 0 0 0
final struc =
rtce= 0 zrce= 0 fnce= 5 ince= 5
rtre= 0 zrre= 0 fnre= 5 inre= 5
computed eigenvalues
eigenvalue= 0.12317D+02 0.00000D+00
eigenvalue= 0.90161D+01 0.00000D+00
eigenvalue= 0.60000D+01 0.00000D+00
eigenvalue= 0.29839D+01 0.00000D+00
eigenvalue= -0.31688D+00 0.00000D+00
guptri - workspace for rzstr - 2 2 2 2 1
3 7 10 22 12 14 26 16 18 20
guptri - m,n,epsu= 2 2 0.222045D-13
reduction 1
kstr, last= 1
1
0
0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 2
1 2
-1 0
-1 0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
final kstr=
1 2 3 4
-1 -1 2 -1
-1 -1 2 -1
nsumrz= 0
rsumrz= 0
djordz= 0
nsumli= 0
rsumli= 0
djordi= 0
dimreg= 2
final pstruc= 0 0 0 0
final struc =
rtce= 0 zrce= 0 fnce= 2 ince= 2
rtre= 0 zrre= 0 fnre= 2 inre= 2
computed eigenvalues
eigenvalue= 0.65865D-01 0.00000D+00
eigenvalue= -0.65865D-01 0.00000D+00
guptri - workspace for rzstr - 3 3 3 3 1
4 13 17 35 20 23 44 26 29 32
guptri - m,n,epsu= 3 3 0.222045D-13
reduction 1
kstr, last= 1
1
1
0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 3
1 2 3
1 -1 1
0 -1 0
accumulated perturbations in a,b = 0.634201D-15 0.000000D+00
final kstr=
1 2 3 4 5 6
1 -1 1 -1 2 -1
0 -1 0 -1 2 -1
nsumrz= 1
rsumrz= 0
djordz= 0
nsumli= 1
rsumli= 0
djordi= 0
dimreg= 2
final pstruc= 1 1 1 2
final struc =
1 1
rtce= 1 zrce= 1 fnce= 3 ince= 3
rtre= 0 zrre= 0 fnre= 2 inre= 2
computed eigenvalues
eigenvalue= -0.65865D-01 0.00000D+00
eigenvalue= 0.65865D-01 0.00000D+00
guptri - workspace for rzstr - 4 4 4 4 1
5 21 26 50 30 34 66 38 42 46
guptri - m,n,epsu= 4 4 0.222045D-13
reduction 1
kstr, last= 1
1
2
0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 3
1 2 3
2 -1 2
0 -1 0
accumulated perturbations in a,b = 0.634201D-15 0.000000D+00
final kstr=
1 2 3 4 5 6
2 -1 2 -1 2 -1
0 -1 0 -1 2 -1
nsumrz= 2
rsumrz= 0
djordz= 0
nsumli= 2
rsumli= 0
djordi= 0
dimreg= 2
final pstruc= 1 1 1 2
final struc =
2 2
rtce= 2 zrce= 2 fnce= 4 ince= 4
rtre= 0 zrre= 0 fnre= 2 inre= 2
computed eigenvalues
eigenvalue= -0.65865D-01 0.00000D+00
eigenvalue= 0.65865D-01 0.00000D+00
guptri - workspace for rzstr - 5 5 5 5 1
6 31 37 67 42 47 92 52 57 62
guptri - m,n,epsu= 5 5 0.222045D-13
reduction 1
kstr, last= 1
1
3
0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 3
1 2 3
3 -1 3
0 -1 0
accumulated perturbations in a,b = 0.634201D-15 0.000000D+00
final kstr=
1 2 3 4 5 6
3 -1 3 -1 2 -1
0 -1 0 -1 2 -1
nsumrz= 3
rsumrz= 0
djordz= 0
nsumli= 3
rsumli= 0
guptri_: info=0
guptri_: rtre=0, rtce=3, zrre=0, zrce=3, fnre=2, fnce=5, inre=2, ince=5
guptri_: dim regular part = 2, zero mult = 0, infty mult = 0, total mult = 2
eigenvalues before reordering:
eigenvalue = (-0.0658653, 0)
eigenvalue = (0.0658653, 0)
Mihelcic4 5x5 output
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.0565, -0.0000) ( -3.6535, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0162, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
Mihelcic4 5x5 output
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.8583, 0.0000) ( -0.4566, -0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.2457, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( -0.0000, 0.0000) ( -0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000)
djordi= 0
dimreg= 2
final pstruc= 1 1 1 2
final struc =
3 3
rtce= 3 zrce= 3 fnce= 5 ince= 5
rtre= 0 zrre= 0 fnre= 2 inre= 2
computed eigenvalues
eigenvalue= -0.65865D-01 0.00000D+00
eigenvalue= 0.65865D-01 0.00000D+00
guptri - workspace for rzstr - 2 2 2 2 1
3 7 10 22 12 14 26 16 18 20
guptri - m,n,epsu= 2 2 0.222045D-13
reduction 1
kstr, last= 2
1 2
1 0
1 0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
reduction 4, kfirst= 4
1 2 3 4
1 0 -1 0
1 0 -1 0
accumulated perturbations in a,b = 0.000000D+00 0.000000D+00
final kstr=
1 2 3 4 5 6
1 0 -1 -1 1 -1
1 0 -1 -1 1 -1
nsumrz= 1
rsumrz= 1
djordz= 1
nsumli= 0
rsumli= 0
djordi= 0
dimreg= 1
final pstruc= 0 1 1 1
final struc =
1
rtce= 0 zrce= 1 fnce= 2 ince= 2
rtre= 0 zrre= 1 fnre= 2 inre= 2
computed eigenvalues
eigenvalue= -0.00000D+00 -0.00000D+00
eigenvalue= 0.10000D+01 -0.00000D+00
guptri_: info=0
guptri_: rtre=0, rtce=0, zrre=1, zrce=1, fnre=2, fnce=2, inre=2, ince=2
guptri_: dim regular part = 2, zero mult = 1, infty mult = 0, total mult = 1
eigenvalues before reordering:
eigenvalue = (-0, -0)
eigenvalue = (1, -0)