@@ -2103,122 +2103,6 @@ ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A)
21032103 return eigenpair;
21042104}
21052105
2106- std::vector<ComplexEigenPair> NonSymmetricEigenSolve (const Matrix & A)
2107- {
2108- char jobvl = ' V' , jobrl = ' V' ;
2109-
2110- int info;
2111- int k = A.numColumns ();
2112- int lwork = std::max (k*k, 10 *k);
2113- double * work = new double [lwork];
2114- double * e_real = new double [k];
2115- double * e_imaginary = new double [k];
2116- Matrix* ev_l = new Matrix (k, k, false );
2117- Matrix* ev_r = new Matrix (k, k, false );
2118- Matrix* A_copy = new Matrix (A);
2119-
2120- // A now in a row major representation. Put it
2121- // into column major order.
2122- for (int row = 0 ; row < k; ++row) {
2123- for (int col = row+1 ; col < k; ++col) {
2124- double tmp = A_copy->item (row, col);
2125- A_copy->item (row, col) = A_copy->item (col, row);
2126- A_copy->item (col, row) = tmp;
2127- }
2128- }
2129-
2130- // Now call lapack to do the eigensolve.
2131- dgeev (&jobvl, &jobrl, &k, A_copy->getData (), &k, e_real, e_imaginary, ev_l->getData (),
2132- &k, ev_r->getData (), &k, work, &lwork, &info);
2133-
2134- // Eigenvectors now in a column major representation. Put it
2135- // into row major order.
2136- for (int row = 0 ; row < k; ++row) {
2137- for (int col = row+1 ; col < k; ++col) {
2138- double tmp = ev_r->item (row, col);
2139- ev_r->item (row, col) = ev_r->item (col, row);
2140- ev_r->item (col, row) = tmp;
2141-
2142- tmp = ev_l->item (row, col);
2143- ev_l->item (row, col) = ev_l->item (col,row);
2144- ev_l->item (col,row) = tmp;
2145- }
2146- }
2147-
2148- std::vector<ComplexEigenPair> eigenpairs;
2149-
2150- ComplexEigenPair eigenpair_right;
2151- eigenpair_right.ev_real = new Matrix (k, k, false );
2152- eigenpair_right.ev_imaginary = new Matrix (k, k, false );
2153-
2154- // Separate lapack eigenvector into real and imaginary parts
2155- for (int i = 0 ; i < k; ++i)
2156- {
2157- for (int row = 0 ; row < k; ++row) {
2158- eigenpair_right.ev_real ->item (row, i) = ev_r->item (row, i);
2159- }
2160- if (e_imaginary[i] != 0 )
2161- {
2162- for (int row = 0 ; row < k; ++row) {
2163- eigenpair_right.ev_real ->item (row, i + 1 ) = ev_r->item (row, i);
2164- eigenpair_right.ev_imaginary ->item (row, i) = ev_r->item (row, i + 1 );
2165- eigenpair_right.ev_imaginary ->item (row, i + 1 ) = -ev_r->item (row, i + 1 );
2166- }
2167-
2168- // Skip the next eigenvalue since it'll be part of the complex
2169- // conjugate pair.
2170- ++i;
2171- }
2172- }
2173-
2174- for (int i = 0 ; i < k; i++)
2175- {
2176- eigenpair_right.eigs .push_back (std::complex <double >(e_real[i], e_imaginary[i]));
2177- }
2178-
2179- eigenpairs.push_back (eigenpair_right);
2180-
2181- ComplexEigenPair eigenpair_left;
2182- eigenpair_left.ev_real = new Matrix (k, k, false );
2183- eigenpair_left.ev_imaginary = new Matrix (k, k, false );
2184-
2185- // Separate lapack eigenvector into real and imaginary parts
2186- for (int i = 0 ; i < k; ++i)
2187- {
2188- for (int row = 0 ; row < k; ++row) {
2189- eigenpair_left.ev_real ->item (row, i) = ev_l->item (row, i);
2190- }
2191- if (e_imaginary[i] != 0 )
2192- {
2193- for (int row = 0 ; row < k; ++row) {
2194- eigenpair_left.ev_real ->item (row, i + 1 ) = ev_l->item (row, i);
2195- eigenpair_left.ev_imaginary ->item (row, i) = ev_l->item (row, i + 1 );
2196- eigenpair_left.ev_imaginary ->item (row, i + 1 ) = -ev_l->item (row, i + 1 );
2197- }
2198-
2199- // Skip the next eigenvalue since it'll be part of the complex
2200- // conjugate pair.
2201- ++i;
2202- }
2203- }
2204-
2205- for (int i = 0 ; i < k; i++)
2206- {
2207- eigenpair_left.eigs .push_back (std::complex <double >(e_real[i], e_imaginary[i]));
2208- }
2209-
2210- eigenpairs.push_back (eigenpair_left);
2211-
2212- delete [] work;
2213- delete [] e_real;
2214- delete [] e_imaginary;
2215- delete ev_r;
2216- delete ev_l;
2217- delete A_copy;
2218-
2219- return eigenpairs;
2220- }
2221-
22222106void SerialSVD (Matrix* A,
22232107 Matrix* U,
22242108 Vector* S,
0 commit comments