DMelt:Numeric/4 Linear Equations

From HandWiki
Member

Linear Equations

Linear_equation is an algebraic equation in which each term is either a constant or the product of a constant and (the first power of) a single variable. A linear system is any system than can be expressed in the format Ax = b. where A is m by n, x is n by o, and b is m by o. Most of the time o=1. There are numerous ways to solve such system using DMelt.

Solving linear systems

First we consider the Apache math library. Consider a linear systems of equations of the form AX=B. For example, consider 

2x + 3y - 2z =  1
   x + 7y + 6x = -2
  4x - 3y - 5z =  1

We will solve this using DecompositionSolver DecompositionSolver of the Apache Common Math package: 

from org.apache.commons.math3.linear import *

# get the coefficient matrix A using LU decomposition
coeff= Array2DRowRealMatrix([[2,3,-2],[-1,7,6],[4,-3,-5]])
solver =LUDecomposition(coeff).getSolver()

#  use solve(RealVector) to solve the system 
constants = ArrayRealVector([1, -2, 1 ])
solution = solver.solve(constants);

print "Solution: x=",solution.getEntry(0), "y=",solution.getEntry(1),"z=",solution.getEntry(2)

The execution of this code prints: 

Solution: x= -0.369863013699 y= 0.178082191781 z= -0.602739726027

Read more for different types of decomposition here.

Using EJML library

A generic way to solve linear equations using EJML is shown below: 

A = DenseMatrix64F(m,n);
x = DenseMatrix64F(n,1);
b = DenseMatrix64F(m,1);

.... code to fill matrices ....

if !CommonOps.solve(A,b,x) : 
     print "Singular matrix";

Linear solvers will in general fail (with some notable exceptions) to produce a meaning full solution if the 'A' matrix is singular. When 'A' is singular then there is an infinite number of solutions. Below is a Jython code which solves a system of linear equations: 

from org.ejml.data import DenseMatrix64F 
from org.ejml.ops import RandomMatrices 
from org.ejml.factory import LinearSolverFactory 
from java.util import Random
rand =Random()

A=DenseMatrix64F(3,3, True, [5, 2, 3, 1.5, -2, 8, -3, 4.7, -0.5])
b=DenseMatrix64F(3,1, True, [18, 21.5, 4.9000])
x=RandomMatrices.createRandom(3,1,rand)
solver =LinearSolverFactory.linear(3)

if   solver.setA(A) == False:
   print "Singular matrix"

c=solver.solve(b,x)
x.print()

Solving equation using multithreaded approach

You can solve linear equations using multiple cores. This addresses the needs of high performance computing (HPC) community. In this approach, the calculations are significantly faster (roughly proportional to the number of available cores of the computer).


  • Solving a system of linear equations using QR factorization.
  • Solving a system of linear equations using LU factorization.
  • Solving a system of linear equations using Cholesky factorization.

Here is the example for solving a system of linear equations using LU factorization using 4 threads (cores). If you do not specify the numbers of cores, it the program will use the maximal number of available cores seen by JVM.

Example for solving a system of linear equations using LU factorization using multiple cores. We will use Dplasma Dplasma class. 

from  edu.emory.mathcs.jplasma.tdouble import Dplasma;
from  edu.emory.mathcs.utils import ConcurrencyUtils;
from jarray import *
from java.lang import Math

M=15;	N=10;	NRHS=5;
A=zeros(M*N, 'd'); 	B=zeros(M*NRHS,'d')

for i in range(M):
	for j in range(N):
		A[M*j + i] = 0.5 - Math.random()

for i  in range(M):
	for  j  in range (NRHS):
		B[M * j + i] = Math.random()

ConcurrencyUtils.setNumberOfThreads(2) # (Optional step) Set the maximum number of threads
Dplasma.plasma_Init(M, N, NRHS)              #  Plasma Initialize
Dplasma.plasma_set_int(Dplasma.PLASMA_CONCURRENCY, 2)

L=Dplasma.plasma_Allocate_L(N, N);
IPIV=Dplasma.plasma_Allocate_IPIV(N, N);
INFO = Dplasma.plasma_DGESV(N, NRHS, A, 0, N, L, 0, IPIV, 0, B, 0, N);
Dplasma.plasma_Finalize()
if INFO < 0:
            print "-- Error in DgelsExample example !"
else:
            print  "-- Run successfull !"
            print IPIV

This is example for the Cholesky method on multiple cores: 

from  edu.emory.mathcs.jplasma.tdouble import Dplasma;
from  edu.emory.mathcs.utils import ConcurrencyUtils;
from jarray import *
from java.lang import Math

M=15;	N=10;	NRHS=5;
A=zeros(M*N, 'd'); 	B=zeros(M*NRHS,'d')

for i in range(M):
	for j in range(N):
		A[M*j + i] = 0.5 - Math.random()

for i  in range(M):
	for  j  in range (NRHS):
		B[M * j + i] = Math.random()

ConcurrencyUtils.setNumberOfThreads(4) # (Optional step) Set the maximum number of threads
Dplasma.plasma_Init(M, N, NRHS)               #  Plasma Initialize
Dplasma.plasma_set_int(Dplasma.PLASMA_CONCURRENCY, 2)
INFO = Dplasma.plasma_DPOSV(Dplasma.PlasmaUpper, N, NRHS, A, 0, N, B, 0, N);
Dplasma.plasma_Finalize()
if INFO < 0:
            print "-- Error in DgelsExample example !"
else:
            print  "-- Run successful !"
Retrieved from "https://handwiki.org/wiki/index.php?title=DMelt:Numeric/4_Linear_Equations&oldid=220"