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


DevX Home    Today's Headlines   Articles Archive   Tip Bank   Forums   

Results 1 to 3 of 3

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

Hybrid View

  1. #1
    parag Guest

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


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

  2. #2
    jonnin Guest

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


    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++.



  3. #3
    Master Flex Guest

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


    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)

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  
HTML5 Development Center
 
 
FAQ
Latest Articles
Java
.NET
XML
Database
Enterprise
Questions? Contact us.
C++
Web Development
Wireless
Latest Tips
Open Source


   Development Centers

   -- Android Development Center
   -- Cloud Development Project Center
   -- HTML5 Development Center
   -- Windows Mobile Development Center