Source Code for Module minimise.hessian_mods.gmw81_old

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003 Edward d'Auvergne                                        # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23   
 24  from LinearAlgebra import eigenvectors, inverse 
 25  from Numeric import Float64, array, dot, matrixmultiply, sort, sqrt, transpose 
 26   
 27   
 28 -def gmw_old(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0): 
 29      """Modified Cholesky Hessian modification. 
 30   
 31      Algorithm 6.5 from page 148. 
 32   
 33      Somehow the data structures l, d, and c can be stored in d2fk! 
 34      """ 
 35   
 36      # Debugging (REMOVE!!!). 
 37      #d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], Float64) 
 38   
 39      if print_flag >= 3: 
 40          print "d2fk: " + `d2fk` 
 41          eigen = eigenvectors(d2fk) 
 42          eigenvals = sort(eigen[0]) 
 43          print "Eigenvalues: " + `eigenvals` 
 44   
 45      # Calculate gamma(A) and xi(A). 
 46      gamma = 0.0 
 47      xi = 0.0 
 48      for i in xrange(n): 
 49          gamma = max(abs(d2fk[i, i]), gamma) 
 50          for j in xrange(i+1, n): 
 51              xi = max(abs(d2fk[i, j]), xi) 
 52   
 53      # Calculate delta and beta. 
 54      delta = mach_acc * max(gamma + xi, 1) 
 55      beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc)) 
 56   
 57      # Initialise data structures. 
 58      d2fk_orig = 1.0 * d2fk 
 59      c = 0.0 * d2fk 
 60      d = 0.0 * dfk 
 61      l = 1.0 * I 
 62      e = 0.0 * dfk 
 63      P = 1.0 * I 
 64   
 65      # Initial diagonal elements of c. 
 66      for k in xrange(n): 
 67          c[k, k] = d2fk[k, k] 
 68   
 69      # Main loop. 
 70      for j in xrange(n): 
 71          if print_flag >= 3: 
 72              print "\n<j: " + `j` + ">" 
 73   
 74          # Row and column swapping. 
 75          p = 1.0 * I 
 76          q = j 
 77          for i in xrange(j, n): 
 78              if abs(c[q, q]) <= abs(c[i, i]): 
 79                  q = i 
 80          if print_flag >= 3: 
 81              print "Row and column swapping." 
 82              print "i range: " + `range(j, n)` 
 83              print "q: " + `q` 
 84              print "a: " + `d2fk` 
 85              print "c: " + `c` 
 86              print "d: " + `d` 
 87              print "l: " + `l` 
 88          if q != j: 
 89              # Modify the permutation matrices. 
 90              temp_p, temp_P = 1.0*p[:, q], 1.0*P[:, q] 
 91              p[:, q], P[:, q] = p[:, j], P[:, j] 
 92              p[:, j], P[:, j] = temp_p, temp_P 
 93   
 94              # Permute both c and d2fk. 
 95              c = dot(p, dot(c, p)) 
 96              d2fk = dot(p, dot(d2fk, p)) 
 97   
 98              if print_flag >= 3: 
 99                  print "p: " + `p` 
100                  print "a(mod): " + `d2fk` 
101                  print "c(mod): " + `c` 
102                  print "d(mod): " + `d` 
103                  print "l(mod): " + `l` 
104   
105          # Calculate the elements of l. 
106          if print_flag >= 3: 
107              print "\nCalculate the elements of l." 
108              print "s range: " + `range(j)` 
109          for s in xrange(j): 
110              if print_flag >= 3: 
111                  print "s: " + `s` 
112              l[j, s] = c[j, s] / d[s] 
113          if print_flag >= 3: 
114              print "l: " + `l` 
115   
116          # Calculate c[i, j]. 
117          if print_flag >= 3: 
118              print "\nCalculate c[i, j]." 
119              print "i range: " + `range(j+1, n)` 
120          for i in xrange(j+1, n): 
121              sum = 0.0 
122              for s in xrange(j): 
123                  if print_flag >= 3: 
124                      print "s range: " + `range(j)` 
125                  sum = sum + l[j, s] * c[i, s] 
126              c[i, j] = d2fk[i, j] - sum 
127          if print_flag >= 3: 
128              print "c: " + `c` 
129   
130          # Calculate d. 
131          if print_flag >= 3: 
132              print "\nCalculate d." 
133          theta_j = 0.0 
134          if j < n-1: 
135              if print_flag >= 3: 
136                  print "j < n, " + `j` + " < " + `n-1` 
137                  print "i range: " + `range(j+1, n)` 
138              for i in xrange(j+1, n): 
139                  theta_j = max(theta_j, abs(c[i, j])) 
140          else: 
141              if print_flag >= 3: 
142                  print "j >= n, " + `j` + " >= " + `n-1` 
143          d[j] = max(abs(c[j, j]), (theta_j/beta)**2, delta) 
144          if print_flag >= 3: 
145              print "theta_j: " + `theta_j` 
146              print "d: " + `d` 
147   
148          # Calculate c[i, i]. 
149          if print_flag >= 3: 
150              print "\nCalculate c[i, i]." 
151          if j < n-1: 
152              if print_flag >= 3: 
153                  print "j < n, " + `j` + " < " + `n-1` 
154                  print "i range: " + `range(j+1, n)` 
155              for i in xrange(j+1, n): 
156                  c[i, i] = c[i, i] - c[i, j]**2 / d[j] 
157          else: 
158              if print_flag >= 3: 
159                  print "j >= n, " + `j` + " >= " + `n-1` 
160          if print_flag >= 3: 
161              print "c: " + `c` 
162   
163          # Calculate e. 
164          e[j] = d[j] - c[j, j] 
165   
166      for i in xrange(n): 
167          d2fk[i, i] = d2fk[i, i] + e[i] 
168      d2fk = dot(P, dot(d2fk, transpose(P))) 
169   
170      # Debugging. 
171      if print_flag >= 3: 
172          print "\nFin:" 
173          print "gamma(A): " + `gamma` 
174          print "xi(A): " + `xi` 
175          print "delta: " + `delta` 
176          print "beta: " + `beta` 
177          print "c: " + `c` 
178          print "d: " + `d` 
179          print "l: " + `l` 
180          temp = 0.0 * I 
181          for i in xrange(len(d)): 
182              temp[i, i] = d[i] 
183          print "m: " + `dot(l, sqrt(temp))` 
184          print "mmT: " + `dot(dot(l, sqrt(temp)), transpose(dot(l, sqrt(temp))))` 
185          print "P: " + `P` 
186          print "e: " + `e` 
187          print "d2fk: " + `d2fk` 
188          eigen = eigenvectors(d2fk) 
189          eigenvals = sort(eigen[0]) 
190          print "Eigenvalues: " + `eigenvals` 
191          import sys 
192          sys.exit() 
193   
194      # Calculate the Newton direction. 
195      if return_matrix: 
196          return -matrixmultiply(inverse(d2fk), dfk), d2fk 
197      else: 
198          return -matrixmultiply(inverse(d2fk), dfk) 
199