Porting Matlab to C

Hi chums, I am trying to turn a matlab script into C, and have a problem with one line:

Code:

X = A\b

Where A is a matrix that m~=n and b is a column matrix of proportions m,1.

The full script I am trying to convert (and have up until that point) is:

Code:

%
% gsolve.m − Solve for imaging system response function
%
% Given a set of pixel values observed for several pixels in several
% images with different exposure times, this function returns the
% imaging system’s response function g as well as the log film irradiance
% values for the observed pixels.
%
% Assumes:
%
% Zmin = 0
% Zmax = 255
%
% Arguments:
%
% Z(i,j) is the pixel values of pixel location number i in image j
% B(j) is the log delta t, or log shutter speed, for image j
% l is lamdba, the constant that determines the amount of smoothness
% w(z) is the weighting function value for pixel value z
%
% Returns:
%
% g(z) is the log exposure corresponding to pixel value z
% lE(i) is the log film irradiance at pixel location i
%

function [g,lE]=gsolve(Z,B,l,w)

n = 256;
A = zeros(size(Z,1)*size(Z,2)+n+1,n+size(Z,1));
b = zeros(size(A,1),1);

%% Include the data−fitting equations

k = 1;
for i=1:size(Z,1)
    for j=1:size(Z,2)
        wij = w(Z(i,j)+1);
        A(k,Z(i,j)+1) = wij; A(k,n+i) = −wij; b(k,1) = wij * B(i,j);
        k=k+1;
    end
end

%% Fix the curve by setting its middle value to 0
A(k,129) = 1;
k=k+1;

%% Include the smoothness equations
for i=1:n−2
    A(k,i)=l*w(i+1); A(k,i+1)=−2*l*w(i+1); A(k,i+2)=l*w(i+1);
    k=k+1;
end

%% Solve the system using SVD
x = A\b;


g = x(1:n);
lE = x(n+1:size(x,1));

It is used to generate the response function, using SVD (single value decomposition), of a digital image process by finding the pixel value difference over a set number of known exposure images.

Matlab help defines '\' as:

Code:

\	Backslash or matrix left division. If A is a square matrix, A\B is
roughly the same as inv(A)*B, except it is computed in a different 
way. If A is an n-by-n matrix and B is a column vector with n 
components, or a matrix with several such columns, then X = A\B is 
the solution to the equation AX = B computed by Gaussian 
elimination. A warning message is displayed if A is badly scaled or 
nearly singular. See the reference page for mldivide for more 
information.

	If A is an m-by-n matrix with m ~= n and B is a column vector 
with m components, or a matrix with several such columns, then X =
 A\B is the solution in the least squares sense to the under- or 
overdetermined system of equations AX = B. The effective rank, k, of
 A is determined from the QR decomposition with pivoting (see 
Algorithm for details). A solution X is computed that has at most k 
nonzero components per column. If k < n, this is usually not the 
same solution as pinv(A)*B, which is the least squares solution with
 the smallest norm

Thanks in advance for any help,
Giant.

p.s. If you can think of a better way of doing SVD with C, don't hesitate to tell me , I have run into a total brick wall here....

C has no matrix functions at all, just operations on individual memory addresses. You can emulate matrix functions (after all, Matlab has to do the same thing as most hardware lacks the capability to do matrix math), which is a bit easier in C++.