how to code any simple matrix(system of linear algebra programming) ?

[parag]

how to code any simple matrix(system of linear algebra programming) ?
Please Just send a sample example of matrix programming in c/c++.

[jonnin]

This is unfinished, code I wrote to solve linear algebra ax + xb = c (lyapunov
eqtn). Most of it works, documents tell you what is what. No ramarks on
style and etc, I have 2 comments on that: its ported to ansi c for embedded
(easy port, this is the c++ one), and its hurry up code, < 3 weeks. Has
most common things, inverse, add, transpose, etc.

There is a main in the next post that calls this stuff and how to use it

all.

It was designed to work with Matlab.
//////////////////////////////////////////////////////////////////////////
/*

mainly routines to solve lyapunov eqtn using kronecker products.
other matrix functions added as needed. Image functions added as well.

control type is mask for variable type (int, etc) to make
this reusable.
cryptic array indexing, sorry about that.
define: matrix = one d, row1row2... format.
Note that control_type a[9] works when routines are called with an
array and that control_type * a can be accessed by a[number]; they are
the same thing to C++.

if an existing matrix is assigned a new value from a
routine, the old matrix is a memory leak if not deallocated
manually. (You can do this by creating a tmp pointer and
assigning it to the matrix, then delete [] tmp after the
operation.)

Answers seem to be acceptable except for the do not use sections.
See header section for details.

This is a hack in progress, only a few days work here, deal with it.
most time spent in invert, for speedups.
*/
////////////////////////////////////////////////////////////
// .headers.


#include <iomanip>
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cassert>

using namespace std;



//typedef float control_type; //matrix type
typedef long double control_type; //matrix type

/////////////////////////////////////////////////////////
//rough elapsed times, global
float timer;
void start_time();
void elapsed_time();
/////////////////////////////////////////////////////////

//matrix routines
int checklegal(control_type *m, int n, int o);
//check diagonal for zeros. true if none found.
control_type *ident(int size);
//returns a size * size identity matrix
control_type *kronecker(control_type *a, int n, int m, control_type *b, int
p, int q);
//returns the kronecker product of a and b
control_type *transpose(control_type *m,int r, int c);
//returns the transpose of the matrix
control_type *multiply(control_type *a,int ar,int ac,control_type *b,int
br,int bc);
//returns the standard matrix product. null if bad sizes are sent
control_type *rswap(control_type *m, int n, int o, int r1, int r2);
//returns matrix with r1 and r2 swapped.
control_type *invert(control_type *m, int n, int o);
//returns the matrix inverse, null if cannot do gaussian.
void display(control_type *m, int n, int o);
//prints the matrix to screen. 18 digits, scientific, configurable
control_type *rref(control_type *m, int n, int o);
//returns full row reduced form(identity|answer)
//I am have not updated rref in ages. Probably slow.
control_type *lyapunov(control_type *a,int ar,int ac,control_type *b,int
br, int bc, control_type *c, int cr, int cc);
//returns the solution to the lyapunov general form by kronecker method
//there seem to be many forms of this. I solved
//XA + BX = C. Matlab uses AX+XB = -C. (these are the same but
//matrix C can cause confusion)
control_type *pivotfix(control_type *m, int n, int o);
//does partial pivoting & must keep zeros off diagonal. see above...
//changes sent matrix (identical but new order. assumes scaled
//matrix was sent)
control_type * add(control_type * a,control_type * b, int n, int o);
control_type * sub(control_type * a,control_type * b, int n, int o);
//return sum & difference of a and b matrices.
control_type * scale(control_type *m, int n, int o);
//divide each row by its largest element
control_type * mul_const(control_type *m, int n, int o,control_type con);
//multiply matrix by a constant
control_type * srad(control_type *m, int n, int o);
//spectral radius (largest abs eigenvalue)
//used to find others in other_evals; uses power method.
control_type * other_evals(control_type *m, int n, int o, control_type guess);
//finds eval close to guess using inverse power method
//this is very sensitive to guess!

//both srad and other_evals can find evectors if needed.

inline control_type diff3x3(control_type *m, const int o, const int n, const
long int dx);


/*
this one is for image processing.
The return should go into a matrix as follows:

0's around the edge, as there is no 3x3 on the edge. (returned here)
The max difference between the center and each element touching it
in the other locations. New but works as edge detector.

I WROTE THESE W/O THINKING FIRST. IMAGES ARE COL MAJOR ORDER,
SO 640*480 IS BACKWARDS FOR THIS PROGRAM. So I am reversing row and
col input arg order here and in calls to clean up the interface.
so call them with 480,640 and all is well.

*/
inline control_type mask3x3(control_type *m, const int o, const int n, const
long int dx);
/*
as above, returns:
sum of multiply each element in a 3*3 const matrix by
section of image.


*/


control_type * submatrix(control_type *m, const int mr, const int mc, const
int newr, const int newc, const int dx);
/*
dx = top left corner.
mr,mc = old rows & cols.
Assumes you know the boundry is legal.
returns the new matrix with m[dx] at 0,0 and newr*newc size.

*/

char compare(const control_type *a, const control_type *b, const int n, const
int o, const float tol);
/*
returns 0 if the abs of the max value in a-b > tol,
else returns 1. a and b have same size, n*o

*/



control_type * distort_horizup(control_type * m, int n, int o, float size);
/*
size is a % > 1.0 returns n*size*o matrix. Interpolates values.


*/






















/////////////////////////////////////////////////////////////
//THE FOLLOWING HAVE ISSUES, DO NOT USE w/o looking first.
//some do not even exist.

control_type * javert(control_type * m, int n, int o);
/* an attempt to take the inverse by jacobi method (solving for ax = I)
hope to expand to pseudoinverse and or get a result for matrices that
do not have inverse. (INCOMPLETE!, the srad was 4.1 and must be < 1
to use this. Left as may use sometime for something else.
*/


control_type * sorvert(control_type * m, int n, int o);
/*
continued: trying to get cols of inverse from sor method, using under
relaxtion as above did not converge.
not done yet, determined not faster than gauss.
*/

control_type * qr(control_type * m, int n, int o);
/*
qr factoring -> evalues approx to feed to other_evals...
not done yet
*/


//////////////////////////////////////////////////////////////////////
//.cpp function bodies

control_type *kronecker(control_type *a, int n, int m, control_type *b, int
p, int q)
//returns a kroneckered with b, worst to read, keep going =)
{
control_type *result = new control_type[n*p*m*q]; //return this
int lc1,lc2,lc3,lc4, ndx; // loop counters and index for above

ndx = 0;
//constructed from truth table from defination of problem
for(lc1 = 0; lc1 < n; lc1++) //ar
for(lc2 = 0; lc2 < p; lc2++)//br
for(lc3 = 0; lc3 < m; lc3++)//ac
for(lc4 = 0; lc4 < q; lc4++)//bc
{
//put it into result: row * size of cols + col offset
result[ndx] = a[lc3+((lc1)*m)]*b[((lc2)*q)+lc4];
ndx++; //move in result array
}
return(result);
}

control_type *ident(int size)
{
//returns size * size identity.
control_type *i = new control_type[size*size];
int a,b,dx;
dx = 0;
for (a = 0; a< size; a++)
for(b=0; b<size; b++)
{
if(a == b)
i[dx] = 1;
else
i[dx] = 0;
dx++;
}
return(i);
}


control_type *transpose(control_type *m,int r, int c)
{//returns transpose. user should know new size as
// mxn -> nxm.

//because my matrices are really one D ...
control_type *t = new control_type[r*c];
int j,k,dx;
dx = 0;
for(j = 0; j< c; j++)
for(k =0; k<r; k++)
{
t[dx] = m[k*c+j];
dx++;
}
return(t);
}

control_type *multiply(control_type *a,int ar,int ac,control_type *b,int
br,int bc)
{//returns null for bad dimensions and ab for good ones.

if(ac != br)//then its not valid.
return(NULL);

control_type *product = new control_type[ar*bc];
int i,j,k,dx;

for(i = 0; i<ar*bc;i++)
product[i] = 0; //set up zeros to use += later.

dx = 0;//index into 1d array...

for(i = 0; i<ar;i++)
for(j = 0; j<bc;j++)
{
for(k = 0; k<ac;k++)
{
product[dx]+= a[i*ac+k]*b[k*bc+j];
}
dx++;
}

return(product);

}


control_type *rswap(control_type *m, int n, int o, int r1, int r2)
{//r1,r2 are from row 0
//very fast now

if(r1 == r2) //fast exit for no exchanges
return(m);

int i;
r1 = r1*o; //avoiding the multiply in loop saved time
r2 = r2*o;
control_type tmp;
for(i = 0; i<o; i++)
{
tmp = m[r1+i];
m[r1+i]= m[r2+i];
m[r2+i]=tmp;
}
return(m);
}

int checklegal(control_type *m, int n, int o)
{//true if no zeros on diagonal.
int i;

for( i = 0; i< n; i++) //n*n
{//if(fabs(m[j+i*o])<1e-14) //tol version
if(m[i+i*o]==0)
return(0);
}
return(1);
}


control_type *pivotfix(control_type *m, int n, int o)
{//faster than ever as now seeks good row (1 in proper place)
//before swapping. (thus avoiding more swaps later)
//probably not as fast as it could be, but any lost time is not noticed here
// compared to inverse.


int i,j,kr;
int er;
long double k;

er = n-1;
if(m[er*o+er] == 0) //force last row to not have zero in last col.
{

for(j = 0; j<n-1; j++)
{
if(m[er+j*o] != 0)
{
kr = j;
j = n; //break loop...
}
}
m = rswap(m,n,o,er,kr);

}

//now we don't allow last row to be swapped to fix other rows.
// this may fubar small arrays but skip partial pivot for them
//if its a problem.


for(i = 0; i<n-1; i++)
{
k = 0;
kr = i;
//quit if 1, the max value after scale, is found as best row
if(fabs(m[i+i*o]) != 1) //if row needs a swap
for(j = i+1; j<n-1; j++) //look for good row
{
if(fabs(m[i+j*o]) == 1) //found good row?
{
m = rswap(m,n,o,i,j); //swap good row
j = n+1; //loop breaker
}


else //if not 1 then find and use max remaining.
{
if(fabs(m[i+j*o]) > k)
{

k = fabs(m[i+j*o]);
kr = j;

}
if(j == (n-2))
{
m = rswap(m,n,o,i,kr);
}
}//else
}//for j

}
return(m);
}

control_type *invert(control_type *m, int n, int o)
{//n,o are = for inverse, but for consistency send both.

int i,j,k,dx;

control_type *id;
control_type *result;
control_type *del;
const int o2 = o*2;//strip out multiplies for speed
///////////////////////////////////////////////////
//append identity

id = ident(n);
control_type *new_m = new control_type[o2*n];
dx = 0;
for(i = 0; i<n; i++)
for(j = 0; j<o; j++)
{
new_m[j+i*o2] = m[dx];
new_m[j+i*o2+o] = id[dx];


dx++;

}
///////////////////////////////////////////////////////////
//scale all rows to have max value of 1
del = new_m;
new_m = scale(new_m,n,o2); //new_m gets a new location
delete [] del; //kill old pointer
//partial pivot to maximize diagonal

new_m = pivotfix(new_m,n,o2);


///////////////////////////////////////////////////////////
//row reduce to zeros below diagonal.
for(k = 0;k<n;k++) //for each row
{
const double tmp = 1/ new_m[k*o2+k]; //get scale


for(j = k; j<o2; j++)
{
new_m[k*o2+j] *= tmp; //scale row
}
for(i = k+1; i < n; i++) //for each row under row k
{
const double tmp2 = new_m[i*o2+k]; //get scale

for(j = k+1; j<o2; j++)
{

if(new_m[k*o2+j] != 0)

new_m[i*o2+j] -= new_m[k*o2+j]*tmp2;

//zeroing small elements can help reduce error
//but must choose cutoff carefully.
/*
if(fabs(new_m[i*o2+j]) < 1.0e-15)
new_m[i*o2+j] = 0;
//*/
}


}
}
////////////////////////////////////////////////////////////////
//if not legal then no inverse (non-singular)
if(checklegal(new_m,n,o2) == 0)
return(NULL); //return new_m here for upper triang.

//if you want the results. remove if statement for always
//return upper triangular.
////////////////////////////////////////////////////////////////
// back sub for identity on left, inverse on right.

for(k= o-1; k>0; k--)
for(i=(k-1);i >= 0; i--)//each row b up
{
const double tmp = new_m[k+i*o2];

new_m[i*o2+k] -= tmp;// *new_m[k*o2+k]; // we know is = 1
for (j = o; j< o2; j++) //each col -value
{
new_m[i*o2+j] -= tmp*new_m[k*o2+j];

/*
if(fabs(new_m[i*o2+j]) < 1.0e-15)
new_m[i*o2+j] = 0;
//*/

}
}



//////////////////////////////////////////////////////////////////
//split off and return right side = inverse
result = new control_type[n*o];
for(i = 0; i<n; i++)
for(j = 0; j<n; j++)
result[i*o+j] = new_m[i*o2+j+o];



//clean memory
delete [] id;
delete [] new_m;
return(result);
}


control_type *rref(control_type *m, int n, int o)
{//taken from invert...
// could be improved for speed.

int i,j,k;
control_type tmp;
control_type tmp2;
int o2;
o2 = o;//from above, saves typing.
control_type * new_m;

new_m = scale(m,n,o2);
new_m = pivotfix(new_m,n,o2);

for(k = 0;k<n;k++) //for each row
{
tmp = new_m[k*o2+k]; //get scale
tmp = 1/tmp;

for(j = 0; j<o2; j++)
new_m[k*o2+j] *= tmp; //scale row

for(i = k+1; i < n; i++) //for each row under row k
{
tmp2 = new_m[i*o2+k]; //get scale
for(j = 0; j<o2; j++)
{
if(new_m[k*o2+j] != 0)
new_m[i*o2+j] -= new_m[k*o2+j]*tmp2;

//zeroing small elements can help reduce error
//but must choose cutoff carefully.
/*
if(fabs(new_m[i*o2+j]) < 1.0e-15)
new_m[i*o2+j] = 0;
//*/
}
//subtract scale * row k (scaled to have 1 in first position)

}
}


if(checklegal(new_m,n,o2) == 0)
return(NULL);

for(k= o-1; k>0; k--)
for(i=(k-1);i >= 0; i--)//each row b up
{
tmp = new_m[k+i*o2];
for (j = 0; j< o2; j++) //each col -value
{
if(new_m[k*o2+j] != 0)
new_m[i*o2+j] -= tmp*new_m[k*o2+j];
/*
if(fabs(new_m[i*o2+j]) < 1.0e-15)
new_m[i*o2+j] = 0;
//*/
}
}

return(new_m);
}

control_type * add(control_type * a,control_type * b, int n, int o)
{//adds 2 matrices
int i;
control_type *sum = new control_type[n*o];
for( i = 0; i< n*o; i++)
sum[i] = a[i] + b[i];
return(sum);

}

control_type * sub(control_type * a,control_type * b, int n, int o)
{//subtract 2 matrices

int i;
control_type *sum = new control_type[n*o];
for( i = 0; i< n*o; i++)
sum[i] = a[i] - b[i];
return(sum);

}

void display(control_type *m, int n, int o)
{//display a matrix, "null" if null
//2 options for precision and another for 1's and 0's where
// 1 = valid numbers and 0 = 0's. used to study patterns.

int c;
char ch[50];


cout.setf(ios::scientific);
setprecision(19);

ch[1] = 0;

if(m == NULL)
cout<<"null"<<endl;
else

for(c = 0; c<n*o; c++)
{
sprintf(ch, "%1.19e",m[c]); //compete w/ matlab
//sig figs for humans
// sprintf(ch, "%2.2f",m[c]);
/*

if(m[c] == 0)
ch[0] = '0';
else
ch[0] = '1';
//*/
cout << ch <<" ";
if((c+1) % (o) == 0)
cout<<endl;
}
cout << endl;
}



control_type *lyapunov(control_type *a,int ar,int ac,control_type *b,int
br, int bc, control_type *c, int cr, int cc)
{//lyap by kronecker..

control_type * x;
control_type * cstacked;
control_type * xstacked;

control_type * at = transpose(a,ar,ac);
control_type * in= ident(br);
control_type * im = ident(ar);

control_type * k1;
control_type * k2;

cstacked = transpose(c,cr,cc); //mine already are one d. transpose
//will make matrix stored by cols.


k1 = kronecker(at,ar,ac,in,br,br);
k2 = kronecker(im,ar,ar,b,br,bc);


delete [] im;
delete [] in;

im = add(k1,k2,ar*bc,ar*bc);


in = invert(im,ar*bc,ar*bc);


xstacked = multiply(in,ar*bc,ar*bc,cstacked,ar*bc,1);

delete [] cstacked;
delete [] im;
delete [] in;
delete [] k1;
delete [] k2;
delete [] at;

x = transpose(xstacked,cc,cr);

delete [] xstacked;

return(x);
//result is good, errors from inverse cause minor inaccuracy...

}



control_type * scale(control_type *m, int n, int o)
{//scale largest element in each row to one.

control_type *scaled = new control_type[n*o];
int i,j;
for(i = 0; i<n*o; i++)
scaled[i]=m[i];


control_type s;

for(i = 0; i<n; i++)
{
s = 0;
for(j = 0; j<o; j++)
{
if(fabs(scaled[j+i*o])>s)
s = fabs(scaled[j+i*o]);
}

if(s != 0)//divide check; if 0 then row = 0 and
//multiply is ok.
s = 1/s; //faster to multiply (maybe. hardware has changed)
for(j = 0; j<o; j++)
{
scaled[j+i*o] *= s;
}


}

return(scaled);
}

control_type * mul_const(control_type *m, int n, int o,control_type con)
{//multiple matrix by a scalar
control_type *cp = new control_type[n*o];

int i;
for( i = 0; i<n*o; i++)
cp[i] = m[i]*con;

return(cp);
}


control_type * srad(control_type *m, int n, int o)
{
//computes spectral radius and its evector
//not sure about vector, but srad is ok and in location 0
//of the result

control_type * v;
control_type * u; //always u and v in math because they look alike
//when professors write on board
control_type * r = new control_type[o];
control_type trace;
int i,j;

v = new control_type[o];
u = new control_type[o];
control_type max;

trace = 0;

for(j = 0; j<o; j++)
{
trace += m[j+j*o];
}//sum of diagonal


for( i = 0; i<o; i++)
u[i] = 1e14; //initial guess. how to get a good one?!
//it seems to help if initial guess is > evals.

for( i = 0; i<750; i++)
//too many can crash on memory access... (225*225, 3000 iterations)

{
v = multiply(m,n,o,u,o,1);

max = 0;
for(j = 0; j<o; j++)
{
if(fabs(v[j]) > max)
max = v[j];
max = 1/max;
u = mul_const(v,1,o,max);
}

}


control_type * rt = new control_type;
*rt = v[o-1];

delete [] r;
delete [] u;
delete [] v;

return(rt);

}

control_type * other_evals(control_type *m, int n, int o,control_type guess)
{
control_type * t2;//alternate t2 and t3 as
//sometimes there can be memory errors from calling
//the routines with one var. Windows problem.
//pointless as 220 9433 <3aba55a0$1@news.devx.com> article retrieved - head and body follows
From: "jonnin" <jonnin@vol.com>
Sender: "jonnin" <jonnin@vol.com>
Reply-To: "jonnin" <jonnin@vol.com>
Subject: Re: how to code any simple matrix(system of linear algebra programming) ?
Newsgroups: c++.general
X-User-Info: 207.244.32.59 207.244.32.59
References: <3aba3d64@news.devx.com>
NNTP-Posting-Host: 209.1.14.192
Message-ID: <3aba55a0$1@news.devx.com>
Date: 22 Mar 2001 11:42:24 -0800
X-Trace: 22 Mar 2001 11:42:24 -0800, 209.1.14.192
Lines: 193
Path: news.devx.com
Xref: news.devx.com c++.general:9433


Here is the main to solve an exmaple lyapunov eqtn.
Ugly, but hey, it works...

//////////////////////////////

#include"lyap.cpp"

void main()
{
//////////////////////////////////////////////////////////
//load up matrices to test.

control_type * test;
control_type * test2;

control_type b[9]=
{
10, 1, -1,
1, 15, 1,
-1, 1, 20
};


control_type a[225]={

-1.5950561e+000, 1.0361583e-004, 5.9778891e-002, 1.0082246e+000, -3.8420257e-007,
-3.1718143e-004, 1.1734223e-005, 3.0776635e+001, 4.7353154e-005, 7.9702753e-005,
4.3256539e-004, -1.7269303e-001, 1.5928253e+000, 6.6444270e-001, -5.3041760e-002,
-1.2088862e-002, -2.2753991e+000, 2.2198078e-001, 1.4793196e-004, 1.0115227e+000,
-1.1195844e-003, -3.0813295e+001, 1.0925111e-001, 9.3007793e+001, 2.3720169e-002,
-1.7670620e-007, 5.3531781e-001, 1.4064449e-001, -1.0261129e-001, 1.6838381e-001,
-1.5103893e-001, -2.7868596e-001, -3.4855228e-001, -7.7835957e-002, 1.4013417e-003,
1.0048455e+000, -4.2761221e-002, -9.6409432e+001, 1.2842083e-001, -7.9683372e-005,
3.3545272e-004, -1.2951455e+000, 1.1676357e+001, 4.8092817e+000, -4.5776849e-001,
1.0642631e-001, -1.6393403e-005, -9.7939457e-003, -3.6389801e-004, 4.8009061e-008,
3.3366092e-005, -1.3789586e-006, 3.2105584e+001, -1.6202427e-002, -1.5664713e-002,
-3.0784684e+001, 5.5037894e-002, 6.0775827e-001, 1.4655372e-001, -7.0819329e-002,
2.2981811e-007, 5.6577468e-003, -4.2087370e-006, -1.4270289e-009, -1.9195043e-005,
1.0867894e-008, -3.2084124e+001, -1.0590627e-006, 2.5042339e+000, -1.0898842e+002,
-2.1866742e-010, 3.4092947e+000, -1.0359165e+000, -1.1864058e+000, -2.6954687e+000,
-1.1579689e-010, 2.3782141e-011, -8.2469210e-012, 2.1976768e-011, 1.7250586e-011,
9.4022894e-012, 4.0547044e-009, -2.6883336e-010, 9.2478074e-008, 1.8999998e-002,
9.3474260e+001, -1.7766681e+000, 1.3951397e+001, 6.0867495e+000, -3.1621959e-001,
1.6964488e-006, 2.8153842e-002, -3.1070468e-005, -1.3973072e-008, -1.4166347e-004,
1.0573654e-007, -1.4502086e-004, -1.0316464e-005, -3.3964136e-002, -3.2750502e-001,
-1.1069089e-011, -1.1292761e-002, -7.3484812e-005, 1.7928894e-002, -5.8372296e-002,
-6.5883118e-001, 8.4656958e-005, 6.0571626e-002, 3.3410376e-003, -2.8700259e-007,
-3.0619743e-004, 8.7457689e-006, -9.8429833e-002, -2.9267827e-005, 2.1078812e-007,
1.0011572e+000, -1.0081758e-004, 7.2002273e-002, 3.8535603e-002, 8.9696509e-003,
-5.3641624e-006, -8.9328754e-002, 9.8244378e-005, 4.4106911e-008, 4.4793881e-004,
-3.3376147e-007, -1.2606757e-002, 3.2564277e-005, -2.0135379e+000, 1.0392081e+000,
-4.1645491e-011, -8.6631468e-003, 6.8429668e-004, 1.2008470e-002, 1.1998031e-002,
-1.7177966e-011, -4.2772169e-007, 3.1777834e-010, 1.1727589e-013, 1.4492782e-009,
-8.3883198e-013, -5.9546358e-006, 8.4603729e-011, -9.6649873e-004, 4.9942968e-006,
9.4446259e-012, 1.6564522e-002, -1.0736350e-003, 3.4962406e-003, 2.8923129e-002,
-7.1208234e-013, -5.9718181e-014, -6.8212042e-014, 1.6677557e-013, -1.4918552e-014,
3.7660916e-014, 1.3171983e-011, -4.0918431e-012, 5.6378134e-013, -8.0254205e-012,
-1.0990726e-011, 4.8647835e-003, 2.4219511e-002, 1.8314145e-002, 1.0027723e-002,
-8.8242057e-006, -1.4668658e-001, 1.6161510e-004, 7.2622045e-008, 7.3687269e-004,
-5.4954002e-007, -1.5416376e-001, 5.3617288e-005, -1.5463298e+000, 8.1487962e-003,
-5.4444084e-011, -3.7450586e-001, 2.9513966e-002, 7.2081418e-004, -1.1645637e-001,
-7.0556715e-001, 1.2466662e-005, 8.9414011e-003, 3.5573184e-003, -4.2127594e-008,
-4.4916466e-005, 1.2866208e-006, 1.1283286e-001, -3.6584730e-006, 1.9206809e-007,
7.7130549e-007, 2.9513966e-002, -1.2084138e+000, -4.8825831e-001, 5.4894012e-003,
2.4205254e-002, -4.7177516e-005, -3.3737950e-002, -1.3467437e-004, 1.6002820e-007,
1.7073223e-004, -4.8873551e-006, 8.1289393e-002, 1.4636030e-005, -1.6181713e-008,
-5.1174035e-004, 7.2081419e-004, -4.8825831e-001, -3.3585358e-001, -5.4246275e-002,
-1.2351898e-009, -3.0391532e-005, 2.2615246e-008, 7.0160344e-012, 1.0314210e-007,
-5.7978890e-011, 2.1453499e-001, 5.6642320e-009, 3.2091013e-005, 3.5485157e-004,
-3.4712755e-011, -1.1645637e-001, 5.4894013e-003, -5.4246275e-002, -2.8768219e-001,

};

//for lyap using a
control_type c[225]=
{
4.4566734e+001, 4.0699614e-001, 3.2811884e+002, 1.5485706e+001, -4.3334282e+001,
3.9585986e+002, -3.8328170e-002, 1.7808599e+000, -1.5220553e-001, -3.6264311e-001,
4.3326832e-001, 4.6945062e+000, -2.9395328e+001, -1.1361915e+001, 3.1290102e+000,
4.0699614e-001, 4.9639028e+000, 3.3598321e+000, 1.1036260e+000, 2.0329398e+001,
2.2022936e+000, -3.6335190e-001, -8.7375856e-003, -1.6396274e-001, 8.5968976e-002,
-2.6576934e-002, -1.7244805e+000, -1.4714785e-001, 4.5543733e-001, -3.2733345e-001,
3.2811884e+002, 3.3598321e+000, 2.4158899e+003, 1.1412067e+002, -3.1672447e+002,
2.9144069e+003, -2.9542974e-001, 1.3109341e+001, -1.1342640e+000, -2.6673247e+000,
3.1872158e+000, 3.4425814e+001, -2.1644360e+002, -8.3615716e+001, 2.3042150e+001,
1.5485706e+001, 1.1036260e+000, 1.1412067e+002, 5.7698307e+000, -6.2524151e+000,
1.3663927e+002, -5.5795032e-003, 6.0881051e-001, -9.4444545e-002, -1.3111758e-001,
1.4102457e-001, 1.2358585e+000, -1.0325794e+001, -3.8493772e+000, 1.1953684e+000,
-4.3334282e+001, 2.0329398e+001, -3.1672447e+002, -6.2524151e+000, 3.2100340e+002,
-4.0563945e+002, 1.7047288e+000, -1.9580616e+000, -9.2259968e-001, -1.4276933e-001,
-6.3525945e-001, -1.4619047e+001, 2.5895773e+001, 1.2705548e+001, 2.2722250e+000,
3.9585986e+002, 2.2022936e+000, 2.9144069e+003, 1.3663927e+002, -4.0563945e+002,
3.5188224e+003, -4.7863599e-001, 1.5842774e+001, -1.2746384e+000, -3.1770806e+000,
3.8689306e+000, 4.2362198e+001, -2.6071390e+002, -1.0105694e+002, 2.7345303e+001,
-3.8328170e-002, -3.6335190e-001, -2.9542974e-001, -5.5795032e-003, 1.7047288e+000,
-4.7863599e-001, 8.1387425e-002, -1.6920235e-003, 5.5581971e-003, -2.0604182e-002,
4.5194128e-004, 8.0991854e-002, -3.5802090e-002, -4.2531157e-002, 1.3596269e-001,
1.7808599e+000, -8.7375856e-003, 1.3109341e+001, 6.0881051e-001, -1.9580616e+000,
1.5842774e+001, -1.6920235e-003, 7.3424761e-002, -5.0074126e-003, -1.4369756e-002,
1.7560987e-002, 1.9783407e-001, -1.1718020e+000, -4.5661270e-001, 1.2231414e-001,
-1.5220553e-001, -1.6396274e-001, -1.1342640e+000, -9.4444545e-002, -9.2259968e-001,
-1.2746384e+000, 5.5581971e-003, -5.0074126e-003, 8.6468219e-003, 8.0856890e-005,
-2.9807022e-004, 4.7293298e-002, 1.0310545e-001, 2.1829019e-002, -1.2881680e-002,
-3.6264311e-001, 8.5968976e-002, -2.6673247e+000, -1.3111758e-001, -1.4276933e-001,
-3.1770806e+000, -2.0604182e-002, -1.4369756e-002, 8.0856890e-005, 1.0084915e+000,
-3.6708078e-003, -5.8116510e-002, 2.5638021e-001, 1.0552748e-001, -6.2686926e-002,
4.3326832e-001, -2.6576934e-002, 3.1872158e+000, 1.4102457e-001, -6.3525945e-001,
3.8689306e+000, 4.5194128e-004, 1.7560987e-002, -2.9807022e-004, -3.6708078e-003,
1.0044654e+000, 5.7471941e-002, -2.8410482e-001, -1.1373018e-001, 2.9416514e-002,
4.6945062e+000, -1.7244805e+000, 3.4425814e+001, 1.2358585e+000, -1.4619047e+001,
4.2362198e+001, 8.0991854e-002, 1.9783407e-001, 4.7293298e-002, -5.8116510e-002,
5.7471941e-002, 1.6638836e+000, -3.1049652e+000, -1.3855433e+000, 3.6346505e-001,
-2.9395328e+001, -1.4714785e-001, -2.1644360e+002, -1.0325794e+001, 2.5895773e+001,
-2.6071390e+002, -3.5802090e-002, -1.1718020e+000, 1.0310545e-001, 2.5638021e-001,
-2.8410482e-001, -3.1049652e+000, 1.9983293e+001, 7.5145266e+000, -2.1910657e+000,
-1.1361915e+001, 4.5543733e-001, -8.3615716e+001, -3.8493772e+000, 1.2705548e+001,
-1.0105694e+002, -4.2531157e-002, -4.5661270e-001, 2.1829019e-002, 1.0552748e-001,
-1.1373018e-001, -1.3855433e+000, 7.5145266e+000, 3.4628623e+000, -8.6126403e-001,
3.1290102e+000, -3.2733345e-001, 2.3042150e+001, 1.1953684e+000, 2.2722250e+000,
2.7345303e+001, 1.3596269e-001, 1.2231414e-001, -1.2881680e-002, -6.2686926e-002,
2.9416514e-002, 3.6346505e-001, -2.1910657e+000, -8.6126403e-001, 2.5486729e+001,


};

/*
//orig lyap equation from book
a[0] = 0; a[1] = 1;
a[2] = -2; a[3] = -3;

b[0] = 0; b[1] = 1; b[2] = 0;
b[3] = 0; b[4] = 0; b[5] = 1;
b[6] = -18; b[7] = -27; b[8]= -10;

c[0] = 1; c[1] = 0;
c[2] = 2; c[3] = 1;
c[4] = 3; c[5] = 0;
//*/

/* //evals tester
b[0] = 10; b[1] = -8; b[2] = -4;
b[3] = -8; b[4] = 13; b[5] = 4;
b[6] = -4; b[7] = 5; b[8]= 4;
*/



//////////////////////////////////////////////////////////
//test the functions

test2 = transpose(a,15,15);

// int jlk;
start_time();

test = lyapunov(a,15,15,test2,15,15,c,15,15);


elapsed_time();
//////////////////////////////////////////////////////////
//print results

display(test,15,15);

/*
//col maj to row maj matlab to jon
test = transpose(b,3,2);
display(test,2,3);
//row m to col m jon to matlab
test2 = transpose(test,2,3);
display(test2,1,6);
*/

delete [] test;

}









//////////////////////////////

"parag" <parag_alam@hotmail.com> wrote:
>
>how to code any simple matrix(system of linear algebra programming) ?
>Please Just send a sample example of matrix programming in c/c++.

[Master Flex]

Look for a sample class declaration and class's defined functions and their
use in the book called "Object-Oriented Programming in C++" by Richard Johnsonbaugh,
and Martin Kalin (pp. 301 - 306)