1 // This is core/vnl/algo/vnl_generalized_eigensystem.cxx
2 //
3 // vnl_generalized_eigensystem
4 // Author: Andrew W. Fitzgibbon, Oxford RRG
5 // Created: 29 Aug 96
6 //
7
8 #include <iostream>
9 #include "vnl_generalized_eigensystem.h"
10
11 #include <vnl/vnl_fortran_copy.h>
12 #include <vnl/vnl_matlab_print.h>
13 #include <vnl/algo/vnl_symmetric_eigensystem.h>
14 #include <vnl/algo/vnl_svd.h>
15 #include <vnl/algo/vnl_netlib.h> // rsg_()
16
vnl_generalized_eigensystem(const vnl_matrix<double> & A,const vnl_matrix<double> & B)17 vnl_generalized_eigensystem::vnl_generalized_eigensystem(const vnl_matrix<double>& A,
18 const vnl_matrix<double>& B)
19 :
20 n(A.rows()), V(n,n), D(n)
21 {
22 // Copy source matrices into fortran storage
23 vnl_fortran_copy<double> a(A);
24 vnl_fortran_copy<double> b(B);
25
26 // Make workspace and storage for V transpose
27 vnl_vector<double> work1(n);
28 vnl_vector<double> work2(n);
29 vnl_vector<double> V1(n*n);
30
31 long want_eigenvectors = 1;
32 long ierr = -1;
33
34 // Call EISPACK rsg.
35 v3p_netlib_rsg_ (&n, &n, a, b, D.data_block(),
36 &want_eigenvectors,
37 V1.begin(),
38 work1.begin(),
39 work2.begin(), &ierr);
40
41 // If b was not pos-def, retry with projection onto nullspace
42 if (ierr == 7*n+1) {
43 const double THRESH = 1e-8;
44 vnl_symmetric_eigensystem<double> eig(B);
45 if (eig.D(0,0) < -THRESH) {
46 std::cerr << "**** vnl_generalized_eigensystem: ERROR\n"
47 << "Matrix B is not nonneg-definite\n";
48 vnl_matlab_print(std::cerr, B, "B");
49 std::cerr << "**** eigenvalues(B) = " << eig.D << std::endl;
50 return;
51 }
52 return;
53 }
54
55 // transpose-copy V1 to V
56 {
57 double *vptr = &V1[0];
58 for (int c = 0; c < n; ++c)
59 for (int r = 0; r < n; ++r)
60 V(r,c) = *vptr++;
61 }
62
63 // Diagnose errors
64 if (ierr) {
65 if (ierr == 10*n)
66 std::cerr << "vnl_generalized_eigensystem: N is greater than NM. Bug in interface to rsg.f\n";
67 else {
68 std::cerr << "vnl_generalized_eigensystem: The "
69 << ierr << "-th eigenvalue has not been determined after 30 iterations.\n"
70 << "The eigenvalues should be correct for indices 1.." << ierr-1
71 << ", but no eigenvectors are computed.\n"
72 << "A = " << A
73 << "\nsingular values(A) = " << vnl_svd<double>(A).W() << '\n'
74 << "B = " << B
75 << "\nsingular values(B) = " << vnl_svd<double>(B).W() << '\n';
76 }
77 }
78 }
79