連立一次方程式の解き方 http://qiita.com/fukuroder/items/4b708524783192fc2018

bicgstab.py

# -*- coding: utf-8 -*-
import scipy.io, numpy as np, numpy.linalg as nl

name = 'hcircuit'
eps = 1.0e-3

# read matrix and vector
print '#A:', scipy.io.mminfo(name + '.mtx.gz')
print '#b:', scipy.io.mminfo(name + '_b.mtx.gz')
A = scipy.io.mmread(name + '.mtx.gz').tocsr()
b = scipy.io.mmread(name + '_b.mtx.gz')[:,0]
print '#bicgstab start'

# initialize
x_j = np.zeros(b.size)
r_j = r_0 = p_j = b - A.dot(x_j)
r_0_norm = nl.norm(r_0)

# start
for j in range(1, 100001):
	Ap = A.dot(p_j)
	t = np.dot(r_j,r_0)
	a = t/np.dot(Ap,r_0)
	s = r_j - a*Ap
	As = A.dot(s)
	w = np.dot(As,s)/np.dot(As,As)
	x_j = x_j + a*p_j + w*s
	r_j = s - w*As
	c = np.dot(r_j,r_0)/t*a/w
	p_j = r_j + c*(p_j - w*Ap)

	#convergence test
	#print nl.norm(r_j)/r_0_norm
	if np.dot(r_j,r_j) < eps*eps*r_0_norm*r_0_norm:
		print '#number of iterations:', j
		print '#relative residual: %e' % (nl.norm(r_j)/r_0_norm)
		print '#relative error:    %e' % (nl.norm(b-A.dot(x_j))/r_0_norm)
		break
else:
	print 'do not converge!'