← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/lib/libqt/solve_pep.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libqt/solve_pep.cc
 
1
/*!
 
2
  \file
 
3
  \brief Solve a 2x2 pseudo-eigenvalue problem
 
4
  \ingroup QT
 
5
*/
 
6
 
 
7
#include <cstdio>
 
8
#include <cstdlib>
 
9
#include <cmath>
 
10
 
 
11
extern "C" {
 
12
        
 
13
#define A_MIN 1.0E-10
 
14
 
 
15
/*!
 
16
** solve_2x2_pep(): Solve a 2x2 pseudo-eigenvalue problem of the form
 
17
**    [ H11 - E    H12 - E*S ]  [c1]  
 
18
**    [ H12 - E*S  H22 - E   ]  [c2]  = 0
 
19
**
 
20
**  \param H     =  matrix to get eigenvalues of
 
21
**  \param S     =  overlap between states 1 & 2
 
22
**  \param evals =  pointer to array to hold 2 eigenvalues
 
23
**  \param evecs =  matrix to hold 2 eigenvectors
 
24
**
 
25
** Returns: none
 
26
** \ingroup QT
 
27
*/
 
28
void solve_2x2_pep(double **H, double S, double *evals, double **evecs)
 
29
{
 
30
   int i ;
 
31
   double a, b, c, p, q, x ;
 
32
   double norm, tval, tval2; 
 
33
 
 
34
   /* put in quadratic form */
 
35
   a = 1.0 - S * S ;
 
36
   b = 2.0 * S * H[0][1] - H[0][0] - H[1][1] ;
 
37
   c = H[0][0] * H[1][1] - H[0][1] * H[0][1] ; 
 
38
  
 
39
   /* solve the quadratic equation for E0 and E1 */
 
40
 
 
41
   tval = b*b ;
 
42
   tval -= 4.0 * a * c ;
 
43
   tval2 = sqrt(tval) ;
 
44
   if (tval < 0.0) {
 
45
      fprintf(stderr, "(solve_2x2_pep): radical less than 0\n") ;
 
46
      return ;
 
47
      }
 
48
   else if (fabs(a) < A_MIN) {
 
49
      printf("min a reached\n") ;
 
50
      evals[0] = evals[1] = H[1][1] ;
 
51
      }
 
52
   else { 
 
53
      evals[0] = -b / (2.0 * a) ;
 
54
      evals[0] -= tval2 / (2.0 * a) ; 
 
55
      evals[1] = -b / (2.0 * a) ;
 
56
      evals[1] += tval2 / (2.0 * a) ;
 
57
      }
 
58
 
 
59
   /* Make sure evals[0] < evals[1] */
 
60
   if (evals[1] < evals[0]) { 
 
61
      tval = evals[0] ;  
 
62
      evals[0] = evals[1] ;
 
63
      evals[1] = tval ;
 
64
      } 
 
65
 
 
66
   /* Make sure evals[0] < H[1][1] */
 
67
   if (evals[0] > H[1][1]) {
 
68
      printf("Warning: using H11 as eigenvalues\n") ;
 
69
      evals[0] = evals[1] = H[1][1] ;
 
70
      printf("Got evals[0] = H[1][1] = %12.7f\n", evals[0]) ;
 
71
      }
 
72
 
 
73
   /* get the eigenvectors */
 
74
   for (i=0; i<2; i++) {
 
75
      p = H[0][0] - evals[i] ;
 
76
      q = H[0][1] - S * evals[i] ;
 
77
      x = -p/q ; 
 
78
      norm = 1.0 + x * x ;
 
79
      norm = sqrt(norm) ;
 
80
      evecs[i][0] = 1.0 / norm ;
 
81
      evecs[i][1] = x / norm ;
 
82
      }
 
83
 
 
84
   /* test 
 
85
   for (i=0; i<2; i++) {
 
86
      p = H[i][0] * evecs[0][0] + H[i][1] * evecs[0][1] ;
 
87
      if (i==0) q = evecs[0][0] + S * evecs[0][1] ;
 
88
      else q = S * evecs[0][0] + evecs[0][1] ;
 
89
      q *= evals[0] ;
 
90
      printf("2x2 check %d: LHS = %12.6f  RHS = %12.6f\n", i, p, q) ;
 
91
      }
 
92
    */
 
93
 
 
94
}
 
95
 
 
96
 
 
97
 
 
98
/***
 
99
main()
 
100
{
 
101
double **H, **evecs, *evals ;
 
102
double S ;
 
103
void solve_2x2_pep() ;
 
104
 
 
105
   H = (double **) malloc (2 * sizeof(double *)) ;
 
106
   evecs = (double **) malloc (2 * sizeof(double *)) ;
 
107
   H[0] = (double *) malloc (2 * sizeof(double)) ;
 
108
   H[1] = (double *) malloc (2 * sizeof(double)) ;
 
109
   evecs[0] = (double *) malloc (2 * sizeof(double)) ;
 
110
   evecs[1] = (double *) malloc (2 * sizeof(double)) ;
 
111
   evals = (double *) malloc (2 * sizeof(double)) ;
 
112
   H[0][0] = -83.92663122885 ; H[0][1] = -83.92636151402 ;
 
113
   H[1][0] = -83.92636151402 ; H[1][1] = -83.92663613012 ;
 
114
   S = 0.999996726987 ;
 
115
 
 
116
   solve_2x2_pep(H, S, evals, evecs) ;
 
117
   printf("eval0 = %16.8f\n", evals[0]) ;
 
118
   printf("evec0 = (%12.6f %12.6f)\n", evecs[0][0], evecs[0][1]) ;
 
119
   printf("\neval1 = %16.8f\n", evals[1]) ;
 
120
   printf("evec1 = (%12.6f %12.6f)\n", evecs[1][0], evecs[1][1]) ;
 
121
 
 
122
}
 
123
***/
 
124
 
 
125
} /* extern "C" */