← 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.c

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