Dgefa.py
# -*- coding: utf-8 -*-
"""
Created on Fri May 3 12:09:11 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.Ddot import ddot
#from Documents.ChromaStarPy.GAS.blas.Dscal import dscal
#from Documents.ChromaStarPy.GAS.blas.Idamax import idamax
#from Documents.ChromaStarPy.GAS.blas.Daxpy import daxpy
import Ddot
import Dscal
import Idamax
import Daxpy
def dgefa(a, lda, n):
#a = [ [0.0e0 for i in range(n)] for j in range(n) ] #output array
info = 0
ipvt = [0 for i in range(n)]
#aOut = [ [a[j][i] for i in range(n)] for j in range(n) ]
"""
c
c dgefa factors a double precision matrix by gaussian elimination.
c
c dgefa is usually called by dgeco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c on entry
c
c a double precision(lda, n)
c the matrix to be factored.
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 on return
c
c a an upper triangular matrix and the multipliers
c which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c = 0 normal value.
c = k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that dgesl or dgedi will divide by zero
c if called. use rcond in dgeco for a reliable
c indication of singularity.
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,dscal,idamax
c
c internal variables
c
"""
"""
Port to python by Ian Short
Saint Mary's University
May 2019
"""
#c
#c
#c gaussian elimination with partial pivoting
#c
info = 0
nm1 = n - 1
#print("DGEFA: n ", n, " nm1 ", nm1)
if (nm1 >= 1):
for k in range(nm1):
#print("DGEFA: k ", k, " n-k ", n-k)
kp1 = k + 1
#c
#c find l = pivot index
#c
#l = idamax(n-k+1, a[k][k], 1) + k - 1
#l = idamax(n-k+1, [a[kk][k] for kk in range(k, n)], 1) + k - 1
#print("IDAMAX: a ", [a[kk][k] for kk in range(k, n)])
l = Idamax.idamax(n-k, [a[kk][k] for kk in range(k, n)], 1) + k
#print("l ", l)
ipvt[k] = l
#c
#c zero pivot implies this column already triangularized
#c
#if (a[l][k] != 0.0e0):
if (a[l][k] != 0.0e0):
#c
#c interchange if necessary
#c
if (l != k):
#print("l != k")
#t = a[l][k]
#a[l][k] = a[k][k]
#a[k][k] = t
t = a[l][k]
a[l][k] = a[k][k]
a[k][k] = t
#c
#c compute multipliers
#c
#t = -1.0e0/a[k][k]
t = -1.0e0/a[k][k]
#FORTRAN: call dscal(n-k, t, a[k+1][k], 1)
#3rd parameter is in/out
#a[k+1][k] = dscal(n-k, t, a[k+1][k], 1)
#[a[k+1][kk] for kk in range(k, n)] =\
#dscal(n-k, t, [a[k+1][kk] for kk in range(k, n)], 1)
#print("BEFORE DSCAL: t ", t, " a ", [a[kk][k] for kk in range(k+1, n)])
dscalOut =\
Dscal.dscal(n-k-1, t, [a[kk][k] for kk in range(k+1, n)], 1)
#dscalSize = len(dscalOut)
#[a[k+1][kk] for kk in range(k, n)] = [dscalOut[ll] for ll in range(dscalSize)]
dscalCount = 0
for kk in range(k+1, n):
a[kk][k] = dscalOut[dscalCount]
dscalCount+=1
#print("AFTER DSCAL: a ", [a[kk][k] for kk in range(n)])
#scipy: a[k+1][k] = dscal(t, a[k+1][k], n-k, 1)
#c
#c row elimination with column indexing
#c
for j in range(kp1, n):
#t = a[l][j]
t = a[l][j]
if (l != k):
#a[l][j] = a[k][j]
#a[k][j] = t
a[l][j] = a[k][j]
a[k][j] = t
#FORTRAN call daxpy(n-k, t, a[k+1][k] ,1, a[k+1][j], 1)
#5th parameter is in/out
#a[k+1][j] = daxpy(n-k, t, a[k+1][k] ,1, a[k+1][j], 1)
#[a[k+1][jj] for jj in range(j, n)] =\
#daxpy(n-k, t, [a[k+1][kk] for kk in range(k, n)], 1, [a[k+1][jj] for jj in range(j, n)], 1)
#print("k ", k, " j ", j, " l ", l, " t ", t)
#print("Before DAXPY: [a[kk][j] for kk in range(k+1, n)] ",\
# [a[kk][j] for kk in range(k+1, n)])
daxpyOut =\
Daxpy.daxpy(n-k-1, t, [a[kk][k] for kk in range(k+1, n)], 1, [a[kk][j] for kk in range(k+1, n)], 1)
#daxpySize = len(daxpyOut)
daxpyCount = 0
for kk in range(k+1, n):
a[kk][j] = daxpyOut[daxpyCount]
daxpyCount+=1
#print("After DAXPY: [a[kk][j] for kk in range(k+1, n)] ",\
# [a[kk][j] for kk in range(k+1, n)])
#scipy library: a[k+1][j] = daxpy(t, a[k+1][k], n-k, 1, 1)
if (a[l][k] == 0.0e0):
info = k
#print("DGEFA final n ", n)
#ipvt[n-1] = n
ipvt[n-1] = n-1
if (a[n-1][n-1] == 0.0e0):
#info = n
info = n-1
# Try returning a tupple:
return a, ipvt, info
