ADMB Documentation  Fournier-pthread.1088
 All Classes Namespaces Files Functions Variables Typedefs Friends Defines
eigen_dmat.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: eigen_dmat.cpp 494 2012-06-13 20:41:16Z johnoel $
00003  *
00004  * Copyright (c) 2009-2012 ADMB foundation
00005  */
00011 #include <fvar.hpp>
00012 
00013 void eigens(const dmatrix & m, const dmatrix & _evecs,
00014       const dvector & _evals);
00015 
00023 dvector eigenvalues(const dmatrix &m)
00024 {
00025    if (m.rowsize() != m.colsize())
00026    {
00027       cerr <<
00028    "error -- non square matrix passed to dvector eigenvalues(const dmatrix& m)\n";
00029       ad_exit(1);
00030    }
00031 
00032    int rmin = m.rowmin();
00033    int rmax = m.rowmax();
00034 
00035    dmatrix evecs(rmin, rmax, rmin, rmax);
00036    dvector evals(rmin, rmax);
00037 
00038    eigens(m, evecs, evals);
00039 
00040    return evals;
00041 }
00042 
00050 dmatrix eigenvectors(const dmatrix &m)
00051 {
00052    if (m.rowsize() != m.colsize())
00053    {
00054       cerr <<
00055    "error -- non square matrix passed to dmatrix eigenvectors(const dmatrix& m)\n";
00056       ad_exit(1);
00057    }
00058 
00059    int rmin = m.rowmin();
00060    int rmax = m.rowmax();
00061 
00062    dmatrix evecs(rmin, rmax, rmin, rmax);
00063    dvector evals(rmin, rmax);
00064 
00065    eigens(m, evecs, evals);
00066 
00067    return evecs;
00068 }
00069 
00079 dmatrix eigenvectors(const dmatrix &m, const dvector &_diag)
00080 {
00081    ADUNCONST(dvector, diag) if (m.rowsize() != m.colsize())
00082    {
00083       cerr << "error -- non square matrix passed to dmatrix "
00084    << "eigenvectors(const dmatrix& m, const dvector& _diag)\n";
00085       ad_exit(1);
00086    }
00087 
00088    int rmin = m.rowmin();
00089    int rmax = m.rowmax();
00090 
00091    dmatrix evecs(rmin, rmax, rmin, rmax);
00092    dvector evals(rmin, rmax);
00093    eigens(m, evecs, evals);
00094    diag = evals;
00095 
00096    return evecs;
00097 }
00098 
00099 
00117 void eigens(const dmatrix & m, const dmatrix & _evecs,
00118       const dvector & _evals)
00119 {
00120    ADUNCONST(dmatrix, evecs)
00121       ADUNCONST(dvector, evals) if (m.rowsize() != m.colsize())
00122    {
00123       cerr << "error -- non square matrix passed to "
00124    <<
00125    "void (const dmatirx& m, const dmatirx& evecs, const dvector& evals)\n";
00126       ad_exit(1);
00127    }
00128    if ((m.rowsize() != evecs.rowsize())
00129        || (m.colsize() != evecs.colsize()))
00130    {
00131       cerr << "error -- eigenvector array not compatible with matrix in "
00132    <<
00133    "void (const dmatirx& m, const dmatirx& evecs, const dvector& evals)\n";
00134       ad_exit(1);
00135    }
00136    if (m.rowsize() != evals.size())
00137    {
00138       cerr << "error -- eigenvalue array not compatible with matrix in "
00139    <<
00140    "void (const dmatirx& m, const dmatirx& evecs, const dvector& evals)\n";
00141       ad_exit(1);
00142    }
00143 
00144    int evecslb = evecs.rowmin();
00145    int evalslb = evals.indexmin();
00146 
00147    int N = m.colsize();
00148 
00149    dmatrix mat = symmetrize(m);
00150 
00151    //shift m, evecs, and evals so they are compatible with algorithm
00152    mat.rowshift(0);
00153    mat.colshift(0);
00154    evecs.rowshift(0);
00155    evecs.colshift(0);
00156    evals.shift(0);
00157 
00158 
00159    //mat[ row, column ] = A[ (row*row+row)/2 + column ]
00160    dvector A(0, N * (N + 1) / 2);
00161 
00162    for (int i = 0; i < N; i++)
00163    {
00164       for (int j = 0; j <= i; j++)
00165       {
00166    A((i * i + i) / 2 + j) = mat(i, j);
00167       }
00168    }
00169 
00170    dvector RR(0, N * N - 1);
00171    dvector E(0, N - 1);
00172 
00173    // run cephes code
00174    int IND, L, LL, LM, M, MM, MQ, I, J, IA, LQ;
00175    int IQ, IM, IL, NLI, NMI;
00176    double ANORM, ANORMX, AIA, THR, ALM, ALL, AMM, X, Y;
00177    double SINX, SINX2, COSX, COSX2, SINCS, AIL, AIM;
00178    double RLI, RMI;
00179    static double RANGE = 1.0e-10; /*3.0517578e-5; */
00180 
00181 
00182    /* Initialize identity matrix in RR[] */
00183    for (J = 0; J < N * N; J++)
00184       RR[J] = 0.0;
00185    MM = 0;
00186    for (J = 0; J < N; J++)
00187    {
00188       RR[MM + J] = 1.0;
00189       MM += N;
00190    }
00191 
00192    ANORM = 0.0;
00193    for (I = 0; I < N; I++)
00194    {
00195       for (J = 0; J < N; J++)
00196       {
00197    if (I != J)
00198    {
00199       IA = I + (J * J + J) / 2;
00200       AIA = A[IA];
00201       ANORM += AIA * AIA;
00202    }
00203       }
00204    }
00205    if (ANORM <= 0.0)
00206       goto done;
00207    ANORM = sqrt(ANORM + ANORM);
00208    ANORMX = ANORM * RANGE / N;
00209    THR = ANORM;
00210 
00211    while (THR > ANORMX)
00212    {
00213       THR = THR / N;
00214 
00215       do
00216       {       /* while IND != 0 */
00217    IND = 0;
00218 
00219    for (L = 0; L < N - 1; L++)
00220    {
00221 
00222       for (M = L + 1; M < N; M++)
00223       {
00224          MQ = (M * M + M) / 2;
00225          LM = L + MQ;
00226          ALM = A[LM];
00227          if (fabs(ALM) < THR)
00228       continue;
00229 
00230          IND = 1;
00231          LQ = (L * L + L) / 2;
00232          LL = L + LQ;
00233          MM = M + MQ;
00234          ALL = A[LL];
00235          AMM = A[MM];
00236          X = (ALL - AMM) / 2.0;
00237          Y = -ALM / sqrt(ALM * ALM + X * X);
00238          if (X < 0.0)
00239       Y = -Y;
00240          SINX = Y / sqrt(2.0 * (1.0 + sqrt(1.0 - Y * Y)));
00241          SINX2 = SINX * SINX;
00242          COSX = sqrt(1.0 - SINX2);
00243          COSX2 = COSX * COSX;
00244          SINCS = SINX * COSX;
00245 
00246          /*       ROTATE L AND M COLUMNS */
00247          for (I = 0; I < N; I++)
00248          {
00249       IQ = (I * I + I) / 2;
00250       if ((I != M) && (I != L))
00251       {
00252          if (I > M)
00253       IM = M + IQ;
00254          else
00255       IM = I + MQ;
00256          if (I >= L)
00257       IL = L + IQ;
00258          else
00259       IL = I + LQ;
00260          AIL = A[IL];
00261          AIM = A[IM];
00262          X = AIL * COSX - AIM * SINX;
00263          A[IM] = AIL * SINX + AIM * COSX;
00264          A[IL] = X;
00265       }
00266       NLI = N * L + I;
00267       NMI = N * M + I;
00268       RLI = RR[NLI];
00269       RMI = RR[NMI];
00270       RR[NLI] = RLI * COSX - RMI * SINX;
00271       RR[NMI] = RLI * SINX + RMI * COSX;
00272          }
00273 
00274          X = 2.0 * ALM * SINCS;
00275          A[LL] = ALL * COSX2 + AMM * SINX2 - X;
00276          A[MM] = ALL * SINX2 + AMM * COSX2 + X;
00277          A[LM] = (ALL - AMM) * SINCS + ALM * (COSX2 - SINX2);
00278       }     /* for M=L+1 to N-1 */
00279    }      /* for L=0 to N-2 */
00280 
00281       }
00282       while (IND != 0);
00283 
00284    }        /* while THR > ANORMX */
00285 
00286  done:;
00287 
00288    /* Extract eigenvalues from the reduced matrix */
00289    L = 0;
00290    for (J = 1; J <= N; J++)
00291    {
00292       L = L + J;
00293       E[J - 1] = A[L - 1];
00294    }
00295 
00296    // end of cephes code
00297 
00298    //put data into evecs and evals
00299    //EV[ n*i+j ] = evecs[i][j]
00300    for (int i = 0; i < N; i++)
00301    {
00302       for (int j = 0; j < N; j++)
00303       {
00304    evecs(j, i) = RR[N * i + j];
00305       }
00306    }
00307 
00308    evals = E;
00309 
00310    //shift evecs and evals back
00311    evecs.rowshift(evecslb);
00312    evecs.colshift(evecslb);
00313    evals.shift(evalslb);
00314 
00315 }