matrix.py

Created by ews31415

Created on September 05, 2022

3.77 KB

Matrix without Numpy


# Math Calculations
#================================
from math import *
#================================
# Matrix Calculations without Numpy
# r: rows, c: colums, M, A, B: Matrices,
# s: scalar

# creating a new matrix
def newmatrix(r,c):
  M = []
  while len(M)<r:
    M.append([])
    while len(M[-1])<c:
      M[-1].append(0.0)
  return M

# M[-1] grabs the last element of M
# Keep in mind that matrices in Python start
# with (0,0) and end in (rows-1,cols-1)
# != is not equals

# Identity Matrix
def identity(n):
  M = newmatrix(n,n)
  for i in range(n):
    M[i][i]=1.0
  return M

# Print a Matrix
def mprint(M):
  for r in M:
    print([x+0 for x in r])

# Transpose a Matrix
def transpose(M):
  # get numbers of rows
  r=len(M)
  # get number of columns
  c=len(M[0])
  # create transpose matrix
  MT=newmatrix(c,r)
  for i in range(r):
    for j in range(c):
      MT[j][i]=M[i][j]
  return MT

# scalar multiplication - EWS
def scalar(M,s):
  r=len(M)
  c=len(M[0])
  MT=newmatrix(r,c)
  for i in range(r):
    for j in range(c):
      MT[i][j]=s*M[i][j]
  return MT

# trace of a square matrix - EWS
def mtrace(M):
  r=len(M)
  c=len(M[0])
  MT=newmatrix(r,c)
  if r !=c:
    return "Not a square matrix"
  else:
    t=0
    for i in range(r):
      t+=M[i][i]
    return t

# matrix addition
def madd(A,B):
  # dimension check
  if len(A) != len(B) or len(A[0]) != len(B[0]):
    return "Dimensions of A and B mismatch"
  else:
    M=newmatrix(len(A),len(A[0]))
    for i in range(len(A)):
      for j in range(len(A[0])):
        M[i][j]=A[i][j]+B[i][j]
    return M

# matrix multiplication
def mmult(A,B):
  # dimension check
  if len(A[0]) != len(B):
    return "Columns of A != Rows of B"
  else:
    M=newmatrix(len(A),len(B[0]))
    for i in range(len(A)):
      for j in range(len(B[0])):
        t=0
        for k in range(len(A[0])):
          t +=A[i][k]*B[k][j]
        M[i][j]=t
    return M

# determinant and inverse of 2 x 2 inverse
# General code to come!
def det2(M):
  if len(M)==2 and len(M[0])==2:
    return M[0][0]*M[1][1]-M[0][1]*M[1][0]
  else:
    return "not a 2 x 2 matrix"

def det3(M):
  if len(M)==3 and len(M[0])==3:
    m10=M[1][0]
    m20=M[2][0]
    m11=M[1][1]
    m21=M[2][1]
    m12=M[1][2]
    m22=M[2][2]
    t=M[0][0]*det2([[m11,m12],[m21,m22]])
    t-=M[0][1]*det2([[m10,m12],[m20,m22]])
    t+=M[0][2]*det2([[m10,m11],[m20,m21]])
    return t
  else:
    return "not a 3 x 3 matrix"

def inv2(M):
  if len(M)==2 and len(M[0])==2:
    t=det2(M)
    MT=newmatrix(2,2)
    MT[0][0]=M[0][0]*1/t
    MT[0][1]=-M[1][0]*1/t
    MT[1][0]=-M[0][1]*1/t
    MT[1][1]=M[1][1]*1/t
    return MT
  else:
    "not a 2 x 2 matrix"

def inv3(M):
  if len(M)==3 and len(M[0])==3:
    t=det3(M)
    A=newmatrix(3,3)
    A[0][0]=det2([[M[1][1],M[1][2]],[M[2][1],M[2][2]]])/t
    A[0][1]=det2([[M[0][2],M[0][1]],[M[2][2],M[2][1]]])/t
    A[0][2]=det2([[M[0][1],M[0][2]],[M[1][1],M[1][2]]])/t
    A[1][0]=det2([[M[1][2],M[1][0]],[M[2][2],M[2][0]]])/t
    A[1][1]=det2([[M[0][0],M[0][2]],[M[2][0],M[2][2]]])/t
    A[1][2]=det2([[M[0][2],M[0][0]],[M[1][2],M[1][0]]])/t
    A[2][0]=det2([[M[1][0],M[1][1]],[M[2][0],M[2][1]]])/t
    A[2][1]=det2([[M[0][1],M[0][0]],[M[2][1],M[2][0]]])/t
    A[2][2]=det2([[M[0][0],M[0][1]],[M[1][0],M[1][1]]])/t
    return A
  else:
    "not a 3 x 3 matrix"

# addition 2022-08-11
# row mutiply and add - this does not change
# the determinant of the original matrix
def rowop(M,r1,s,r2):
  c=len(M[0])
  MT=M
  for k in range (c):
    MT[r2][k]=MT[r1][k]*s+MT[r2][k]
  return MT

# upper triangle of the matrix
# row 0 cannot have a zero
def utri(M):
  MT=M
  c=len(M[0])
  for j in range(c-1):
    for k in range(j+1,c):
      rowop(MT,j,-MT[k][j]/MT[j][j],k)
  return MT

# determinant
def det(M):
  MT=utri(M)
  d=1
  for k in range(len(MT[0])):
    d*=MT[k][k]
  return d