1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/napack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.napack.quasi; 9 10 11 import scid.core.fortran; 12 13 import scid.ports.napack.addchg; 14 import scid.ports.napack.stopit; 15 16 17 18 19 // 20 // ________________________________________________________ 21 // | | 22 // |SOLVE F SUB I (X) = 0, I = 1 TO N, USING A QUASI-NEWTON| 23 // | SCHEME (BROYDEN'S GOOD METHOD) | 24 // | | 25 // | INPUT: | 26 // | | 27 // | X --STARTING GUESS | 28 // | | 29 // | H --STARTING GUESS FOR INVERSE JACOBIAN | 30 // | | 31 // | LH --LEADING (ROW) DIMENSION OF ARRAY H | 32 // | | 33 // | N --NUMBER OF EQUATIONS | 34 // | | 35 // | NDIGIT--DESIRED NUMBER CORRECT DIGITS | 36 // | | 37 // | LIMIT --MAXIMUM NUMBER OF ITERATIONS | 38 // | | 39 // | FUNC --NAME OF FUNCTION EVALUATION SUBROUTINE | 40 // | (EXTERNAL IN MAIN PROG.) FUNC(F,X) PUTS| 41 // | IN ARRAY F THE FUNCTION VALUE AT X | 42 // | | 43 // | W -WORK ARRAY (LENGTH AT LEAST 3N) | 44 // | | 45 // | OUTPUT: | 46 // | | 47 // | X --SOLUTION | 48 // | | 49 // | DIF --INPUT FOR SUBROUTINE WHATIS | 50 // | | 51 // | SIZE --INPUT FOR SUBROUTINE WHATIS | 52 // | | 53 // | PACKAGE SUBROUTINES: ADDCHG,STOPIT | 54 // |________________________________________________________| 55 /// 56 void quasi(Real, Func)(Real[] x_, Real[] h_, int lh, int n, out Real dif, 57 out Real size, int ndigit, int limit, Func func, Real[] w_) 58 { 59 int i, j; 60 Real s, t, u; 61 auto x = dimension(x_.ptr, n); 62 auto h = dimension(h_.ptr, lh, n); 63 auto w = dimension(w_.ptr, n, 3); 64 for (i=1; i<=n; i++) 65 { 66 w[i,1] = 0.0; 67 } 68 goto l30; 69 // -------------------- 70 // |*** ADD S TO X ***| 71 // -------------------- 72 l20:addchg(dif, size, x_, w_, n); 73 stopit(dif, size, ndigit, limit); 74 if (dif <= 0.0) return; 75 // ----------------- 76 // |*** FORM U ***| 77 // ----------------- 78 l30:for (i=1; i<=n; i++) 79 { 80 w[i,2] = w[i,1]; 81 w[i,1] = 0.0; 82 } 83 func(w_[2*lh .. 3*lh], x_); 84 for (j=1; j<=n; j++) 85 { 86 t = w[j,3]; 87 for (i=1; i<=n; i++) 88 { 89 w[i,1] = w[i,1] - t*h[i,j]; 90 } 91 } 92 // ----------------------------------- 93 // |*** FORM S DOT S AND S DOT U ***| 94 // ----------------------------------- 95 t = 0.0; 96 u = 0.0; 97 for (i=1; i<=n; i++) 98 { 99 s = w[i,2]; 100 t = t + s*w[i,1]; 101 u = u + s*s; 102 } 103 t = u - t; 104 if (t == 0.0) goto l20; 105 // ------------------ 106 // |*** UPDATE H ***| 107 // ------------------ 108 for (j=1; j<=n; j++) 109 { 110 s = 0.0; 111 for (i=1; i<=n; i++) 112 { 113 s += h[i,j]*w[i,2]; 114 } 115 s /= t; 116 for (i=1; i<=n; i++) 117 { 118 h[i,j] = h[i,j] + s*w[i,1]; 119 } 120 } 121 u /= t; 122 for (i=1; i<=n; i++) 123 { 124 w[i,1] = u*w[i,1]; 125 } 126 goto l20; 127 } 128 129 130 unittest 131 { 132 alias quasi!(float, void delegate(float[], float[])) fquasi; 133 alias quasi!(double, void delegate(double[], double[])) dquasi; 134 alias quasi!(double, void function(double[], double[])) dfquasi; 135 alias quasi!(real, void delegate(real[], real[])) rquasi; 136 } 137 138 139 version(unittest) 140 { 141 import std.math; 142 import scid.core.testing; 143 } 144 145 unittest 146 { 147 void rosenbrock(double[] fval, double[] x) 148 { 149 fval[0] = (1 - x[0]); 150 fval[1] = 10.0 * (x[1] - x[0]^^2); 151 } 152 153 double[] guess = [-10.0, -5.0]; 154 double[] invJacAtGuess = [-1.0, 20.0, 0.0, 0.1]; 155 double[] work = new double[6]; 156 157 double dif, size; 158 quasi(guess, invJacAtGuess, 2, 2, dif, size, 6, 500, &rosenbrock, work); 159 assert (guess == [1.0, 1.0]); 160 }