1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/linpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost Licence 1.0 7 */ 8 module scid.ports.linpack.gtsl; 9 10 11 import std.math: fabs; 12 13 import scid.core.fortran; 14 15 16 17 18 /// 19 void gtsl(Real)(int n, Real* c_, Real* d_, Real* e_, Real* b_, out int info) 20 { 21 auto c = dimension(c_, n); 22 auto d = dimension(d_, n); 23 auto e = dimension(e_, n); 24 auto b = dimension(b_, n); 25 // 26 // dgtsl given a general tridiagonal matrix and a right hand 27 // side will find the solution. 28 // 29 // on entry 30 // 31 // n integer 32 // is the order of the tridiagonal matrix. 33 // 34 // c double precision(n) 35 // is the subdiagonal of the tridiagonal matrix. 36 // c(2) through c(n) should contain the subdiagonal. 37 // on output c is destroyed. 38 // 39 // d double precision(n) 40 // is the diagonal of the tridiagonal matrix. 41 // on output d is destroyed. 42 // 43 // e double precision(n) 44 // is the superdiagonal of the tridiagonal matrix. 45 // e(1) through e(n-1) should contain the superdiagonal. 46 // on output e is destroyed. 47 // 48 // b double precision(n) 49 // is the right hand side vector. 50 // 51 // on return 52 // 53 // b is the solution vector. 54 // 55 // info integer 56 // = 0 normal value. 57 // = k if the k-th element of the diagonal becomes 58 // exactly zero. the subroutine returns when 59 // this is detected. 60 // 61 // linpack. this version dated 08/14/78 . 62 // jack dongarra, argonne national laboratory. 63 // 64 // no externals 65 // fortran dabs 66 // 67 // internal variables 68 // 69 int k=1,kb=1,kp1=1,nm1=1,nm2=1; 70 Real t; 71 // begin block permitting ...exits to 100 72 // 73 info = 0; 74 c[1] = d[1]; 75 nm1 = n - 1; 76 if (nm1 < 1) goto l40; 77 d[1] = e[1]; 78 e[1] = 0.0; 79 e[n] = 0.0; 80 // 81 for (k=1; k<=nm1; k++) { // end: 30 82 kp1 = k + 1; 83 // 84 // find the largest of the two rows 85 // 86 if (fabs(c[kp1]) < fabs(c[k])) goto l10; 87 // 88 // interchange row 89 // 90 t = c[kp1]; 91 c[kp1] = c[k]; 92 c[k] = t; 93 t = d[kp1]; 94 d[kp1] = d[k]; 95 d[k] = t; 96 t = e[kp1]; 97 e[kp1] = e[k]; 98 e[k] = t; 99 t = b[kp1]; 100 b[kp1] = b[k]; 101 b[k] = t; 102 l10: ; 103 // 104 // zero elements 105 // 106 if (c[k] != 0.0) goto l20; 107 info = k; 108 // ............exit 109 goto l100; 110 l20: ; 111 t = -c[kp1]/c[k]; 112 c[kp1] = d[kp1] + t*d[k]; 113 d[kp1] = e[kp1] + t*e[k]; 114 e[kp1] = 0.0; 115 b[kp1] = b[kp1] + t*b[k]; 116 l30: ;} 117 l40: ; 118 if (c[n] != 0.0) goto l50; 119 info = n; 120 goto l90; 121 l50: ; 122 // 123 // back solve 124 // 125 nm2 = n - 2; 126 b[n] = b[n]/c[n]; 127 if (n == 1) goto l80; 128 b[nm1] = (b[nm1] - d[nm1]*b[n])/c[nm1]; 129 if (nm2 < 1) goto l70; 130 for (kb=1; kb<=nm2; kb++) { // end: 60 131 k = nm2 - kb + 1; 132 b[k] = (b[k] - d[k]*b[k+1] - e[k]*b[k+2])/c[k]; 133 l60: ;} 134 l70: ; 135 l80: ; 136 l90: ; 137 l100: ; 138 // 139 return; 140 } 141 142 143 unittest 144 { 145 alias gtsl!float fgtsl; 146 alias gtsl!double dgtsl; 147 alias gtsl!real rgtsl; 148 }