https://github.com/sevenian3/ChromaStarPy
Tip revision: 103d3d0df6d9574c49818f149e1cae8100455d10 authored by Ian Short on 06 July 2023, 18:09:20 UTC
Create ReadMe
Create ReadMe
Tip revision: 103d3d0
Dgesl.py
# -*- coding: utf-8 -*-
"""
Created on Mon May 6 11:42:29 2019
@author:
"""
import math
import numpy
#from scipy.linalg.blas import daxpy
#from scipy.linalg.blas import ddot
#from scipy.linalg.blas import dscal
#from scipy.linalg.blas import idamax
#from Documents.ChromaStarPy.GAS.blas.Daxpy import daxpy
#from Documents.ChromaStarPy.GAS.blas.Ddot import ddot
#from Documents.ChromaStarPy.GAS.blas.Dscal import dscal
#from Documents.ChromaStarPy.GAS.blas.Idamax import idamax
import Daxpy
import Ddot
import Dscal
import Idamax
def dgesl(a, lda, n, ipvt, b, job):
#integer lda,n,ipvt(1),job
#double precision a(lda,1),b(1)
"""
c
c dgesl solves the double precision system
c a * x = b or trans(a) * x = b
c using the factors computed by dgeco or dgefa.
c
c on entry
c
c a double precision(lda, n)
c the output from dgeco or dgefa.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c ipvt integer(n)
c the pivot vector from dgeco or dgefa.
c
c b double precision(n)
c the right hand side vector.
c
c job integer
c = 0 to solve a*x = b ,
c = nonzero to solve trans(a)*x = b where
c trans(a) is the transpose.
c
c on return
c
c b the solution vector x .
c
c error condition
c
c a division by zero will occur if the input factor contains a
c zero on the diagonal. technically this indicates singularity
c but it is often caused by improper arguments or improper
c setting of lda . it will not occur if the subroutines are
c called correctly and if dgeco has set rcond .gt. 0.0
c or dgefa has set info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call dgeco(a,lda,n,ipvt,rcond,z)
c if (rcond is too small) go to ...
c do 10 j = 1, p
c call dgesl(a,lda,n,ipvt,c(1,j),0)
c 10 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,ddot
c
c internal variables
c
"""
#double precision ddot,t
#integer k,kb,l,nm1
#c
nm1 = n - 1
if (job == 0):
#c
#c job = 0 , solve a * x = b
#c first solve l*y = b
#c
if (nm1 >= 1):
for k in range(nm1):
l = ipvt[k]
t = b[l]
if (l != k):
#print("DGESL if triggered")
b[l] = b[k]
b[k] = t
#print("DGESL 1: l ", l, " k, ", k, " b ", b[k])
#FORTRAN call call daxpy(n-k, t, a[k+1][k], 1, b[k+1], 1)
#5th parameter is in/out:
#b[k+1] = daxpy(n-k, t, a[k+1][k], 1, b[k+1], 1)
#[b[kk+1] for kk in range(k, n)] = daxpy(n-k, t,\
# [a[k+1][kk] for kk in range(k, n)], 1, [b[kk+1] for kk in range(k, n)], 1)
daxpyOut =\
Daxpy.daxpy(n-k-1, t, [a[kk][k] for kk in range(k+1, n)], 1, [b[kk] for kk in range(k+1, n)], 1)
daxpyCount = 0
for kk in range(k+1, n):
b[kk] = daxpyOut[daxpyCount]
daxpyCount+=1
#print("DGESL 2: k ", k, " b ", b[k])
#scipy: b[k+1] = daxpy(t, a[k+1][k], n-k, 1, 1)
#c
#c now solve u*x = y
#c
#print("DGESL: Before 2nd DAXPY call n ", n)
for kb in range(n):
#k = n + 1 - kb
k = (n-1) - kb
#print("DGESL: kb ", kb, " k ", k, " b ", b[k], " a ", a[k][k])
b[k] = b[k]/a[k][k]
t = -b[k]
#FORTRAN call: call daxpy(k-1, t, a[1][k], 1, b[1], 1)
#b[1] = daxpy(k-1, t, a[1][k], 1, b[1], 1)
#[b[kk] for kk in range(1, k)] = daxpy(k-1, t,\
# [a[1][kk] for kk in range(1, k)], 1, [b[kk] for kk in range(1, k)], 1)
#print("DGESL: Before DAPXPY 2:")
#print("a ", [a[kk][k] for kk in range(0, k+1)])
#print("b ", [b[kk] for kk in range(0, k+1)])
daxpyOut =\
Daxpy.daxpy(k, t, [a[kk][k] for kk in range(0, k+1)], 1, [b[kk] for kk in range(0, k+1)], 1)
daxpyCount = 0
for kk in range(0, k+1):
b[kk] = daxpyOut[daxpyCount]
daxpyCount+=1
#print("DGESL: After DAPXPY 2:")
#print("b ", [b[kk] for kk in range(0, k+1)])
#scipy: b[0] = daxpy(t, a[0][k], k-1, 1, 1)
# **** goto 100 !!! Oh-oh!!
#c
#c job = nonzero, solve trans(a) * x = b
#c first solve trans(u)*y = b
#c
if (job != 0):
for k in range(n):
#t = ddot(k-1, a[1][k], 1, b[1], 1)
t = Ddot.ddot(k, [a[kk][k] for kk in range(0, k)],\
1, [b[kk] for kk in range(0, k)], 1)
b[k] = (b[k] - t)/a[k][k]
#print("DDOT 1: t ", t)
#c
#c now solve trans(l)*x = y
#c
if (nm1 >= 1):
for kb in range(nm1):
#k = n - kb
k = n - kb - 1
#b[k] = b[k] + ddot(n-k, a[k+1][k], 1, b[k+1], 1)
b[k] = b[k] + Ddot.ddot(n-k, [a[kk][k] for kk in range(k, n)],\
1, [b[kk] for kk in range(k, n)], 1)
#print("DDOT 2: t ", t)
l = ipvt[k]
if (l != k):
t = b[l]
b[l] = b[k]
b[k] = t
return b