standaloneTestFastElim.cpp File Reference

#include <vector>
#include <algorithm>
#include <cassert>
#include <functional>
#include <iostream>
#include <iomanip>
#include "/usr/include/gsl/gsl_linalg.h"
#include "../basecode/SparseMatrix.h"

Include dependency graph for standaloneTestFastElim.cpp:

Go to the source code of this file.

Classes
class	FastElim
class	Unroll

Functions
void	advance (vector< double > &y, const vector< Triplet< double > > &ops, const vector< double > &diagVal)
double	checkAns (const double *m, unsigned int numCompts, const double *ans, const double *rhs)
	main ()
void	sortByColumn (vector< unsigned int > &col, vector< double > &entry)
void	testSorting ()

Variables
const unsigned int	SM_MAX_COLUMNS = 200000
const unsigned int	SM_MAX_ROWS = 200000
const unsigned int	SM_RESERVE = 8

Function Documentation

void advance	(	vector< double > &	y,
		const vector< Triplet< double > > &	ops,
		const vector< double > &	diagVal
	)

Definition at line 402 of file standaloneTestFastElim.cpp.

Referenced by main(), Gsolve::process(), and Ksolve::process().

{
    for ( vector< Triplet< double > >::const_iterator
                i = ops.begin(); i != ops.end(); ++i )
        y[i->c_] -= y[i->b_] * i->a_;

    assert( y.size() == diagVal.size() );
    vector< double >::iterator iy = y.begin();
    for ( vector< double >::const_iterator
                i = diagVal.begin(); i != diagVal.end(); ++i )
        *iy++ *= *i;
}

Here is the caller graph for this function:

double checkAns	(	const double *	m,
		unsigned int	numCompts,
		const double *	ans,
		const double *	rhs
	)

Definition at line 417 of file standaloneTestFastElim.cpp.

Referenced by main().

{
    vector< double > check( numCompts, 0.0 );
    for ( unsigned int i = 0; i < numCompts; ++i ) {
        for ( unsigned int j = 0; j < numCompts; ++j )
            check[i] += m[i*numCompts + j] * ans[j];
    }
    double ret = 0.0;
    for ( unsigned int i = 0; i < numCompts; ++i )
        ret += (check[i] - rhs[i]) * (check[i] - rhs[i] );
    return ret;
}

Here is the caller graph for this function:

main

(

)

Definition at line 434 of file standaloneTestFastElim.cpp.

References advance(), FastElim::buildBackwardSub(), FastElim::buildForwardElim(), checkAns(), SparseMatrix< T >::get(), FastElim::hinesReorder(), FastElim::makeTestMatrix(), and SparseMatrix< T >::print().

{


/*
2    11
 1  4
   3    10
    9  5
     6
     7
     8

        1 2 3 4 5 6 7 8 9 10 11
1       x x x x
2       x x
3       x   x x         x
4       x   x x              x
5               x x     x x
6               x x x   x
7                 x x x
8                   x x
9           x   x x     x
10              x         x
11            x              x
    static double test[] = {
        1,  2,  3,  4,  0,  0,  0,  0,  0,  0,  0,
        5,  6,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        7,  0,  8,  9,  0,  0,  0,  0, 10,  0,  0,
        11, 0, 12, 13,  0,  0,  0,  0,  0,  0, 14,
        0,  0,  0,  0, 15, 16,  0,  0, 17, 18,  0,
        0,  0,  0,  0, 19, 20, 21,  0, 22,  0,  0,
        0,  0,  0,  0,  0, 23, 24, 25,  0,  0,  0,
        0,  0,  0,  0,  0,  0, 26, 27,  0,  0,  0,
        0,  0, 28,  0, 29, 30,  0,  0, 31,  0,  0,
        0,  0,  0,  0, 32,  0,  0,  0,  0, 33,  0,
        0,  0,  0, 34,  0,  0,  0,  0,  0,  0, 35,
    };
    const unsigned int numCompts = 11;
//  static unsigned int parents[] = { 3,1,9,3,6,7,8,-1,6,5,4 };
    static unsigned int parents[] = { 2,0,8,2,5,6,7,-1,5,4,3 };
*/

/*
1   3
 2 4
   7   5
    8 6
     9
     10
     11

        1 2 3 4 5 6 7 8 9 10 11
1       x x
2       x x         x
3           x x
4           x x     x
5               x x
6               x x     x
7         x   x     x x
8                   x x x
9                 x   x x x
10                      x x  x
11                        x  x
    static double test[] = {
        1,  2,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        3,  4,  0,  0,  0,  0,  5,  0,  0,  0,  0,
        0,  0,  6,  7,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  8,  9,  0,  0, 10,  0,  0,  0,  0,
        0,  0,  0,  0, 11, 12,  0,  0,  0,  0,  0,
        0,  0,  0,  0, 13, 14,  0,  0, 15,  0,  0,
        0, 16,  0, 17,  0,  0, 18, 19,  0,  0,  0,
        0,  0,  0,  0,  0,  0, 20, 21, 22,  0,  0,
        0,  0,  0,  0,  0, 23,  0, 24, 25, 26,  0,
        0,  0,  0,  0,  0,  0,  0,  0, 27, 28, 29,
        0,  0,  0,  0,  0,  0,  0,  0,  0, 30, 31,
    };
    const unsigned int numCompts = 11;
    static unsigned int parents[] = { 1,6,3,6,5,8,7,8,9,10,-1};
*/

/*
Linear cable, 12 segments.
*/

    static double test[] = {
        1,  2,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        3,  4,  5,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  6,  7,  8,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  9, 10, 11,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0, 12, 13, 14,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0, 15, 16, 17,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0, 18, 19, 20,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0, 21, 22, 23,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0, 24, 25, 26,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0, 27, 28, 29,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0, 30, 31, 32,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 33, 34,
    };
    const unsigned int numCompts = 12;
    static unsigned int parents[] = { 1,2,3,4,5,6,7,8,9,10,11,-1};

        /*
    static double test[] = {
        1,  2,
        3,  4
    };
    const unsigned int numCompts = 2;
    static double test[] = {
        1, 2, 0, 0,
        3, 4, 5, 0,
        0, 6, 7, 8,
        0, 0, 9, 10
    };
    const unsigned int numCompts = 4;
    static double test[] = {
        1,  2,  0,  0,  0,  0,
        3,  4,  5,  0,  0,  0,
        0,  6,  7,  8,  0,  0,
        0,  0,  9, 10, 11,  0,
        0,  0,  0, 12, 13, 14,
        0,  0,  0,  0, 15, 16,
    };
    const unsigned int numCompts = 6;
    static double test[] = {
        1,  2,  0,  0,  0,  0,
        3,  4,  0,  0,  1,  0,
        0,  0,  7,  8,  0,  0,
        0,  0,  9, 10, 11,  0,
        0,  1,  0, 12, 13, 14,
        0,  0,  0,  0, 15, 16,
    };
    const unsigned int numCompts = 6;
    */
    // testSorting(); // seems to work fine.
    FastElim fe;
    vector< Triplet< double > > fops;
    vector< Unroll > bops;
    fe.makeTestMatrix( test, numCompts );
    fe.print();
    cout << endl << endl;
    vector< unsigned int > parentVoxel;
    parentVoxel.insert( parentVoxel.begin(), &parents[0], &parents[numCompts] );
    fe.hinesReorder( parentVoxel );
    /*
    */
    /*
    vector< unsigned int > shuf;
    for ( unsigned int i = 0; i < numCompts; ++i )
        shuf.push_back( i );
    shuf[0] = 1;
    shuf[1] = 0;
    fe.shuffleRows( shuf );
    */
    fe.print();
    cout << endl << endl;
    FastElim foo = fe;

    vector< unsigned int > diag;
    vector< double > diagVal;
    fe.buildForwardElim( diag, fops );
    fe.print();
    fe.buildBackwardSub( diag, fops, diagVal );
    vector< double > y( numCompts, 1.0 );
    vector< double > ones( numCompts, 1.0 );
    advance( y, fops, diagVal );
    for ( int i = 0; i < numCompts; ++i )
        cout << "y" << i << "]= " << y[i] << endl;

    // Here we verify the answer

    vector< double > alle;
    for( unsigned int i = 0; i < numCompts; ++i ) {
        for( unsigned int j = 0; j < numCompts; ++j ) {
            alle.push_back( foo.get( i, j ) );
        }
    }
    cout << "myCode: " <<
            checkAns( &alle[0], numCompts, &y[0], &ones[0] ) << endl;



    // Here we do the gsl test.
    vector< double > temp( &test[0], &test[numCompts*numCompts] );
    gsl_matrix_view m = gsl_matrix_view_array( &temp[0], numCompts, numCompts );

    vector< double > z( numCompts, 1.0 );
    gsl_vector_view b = gsl_vector_view_array( &z[0], numCompts );
    gsl_vector* x = gsl_vector_alloc( numCompts );
    int s;
    gsl_permutation* p = gsl_permutation_alloc( numCompts );
    gsl_linalg_LU_decomp( &m.matrix, p, &s );
    gsl_linalg_LU_solve( &m.matrix, p, &b.vector, x);
    vector< double > gslAns( numCompts );
    for ( int i = 0; i < numCompts; ++i ) {
        gslAns[i] = gsl_vector_get( x, i );
        cout << "x[" << i << "]=    " << gslAns[i] << endl;
    }
    /*
    */
    cout << "GSL: " << checkAns( test, numCompts, &gslAns[0], &ones[0] ) << endl;
    gsl_vector_free( x );


}

Here is the call graph for this function:

void sortByColumn	(	vector< unsigned int > &	col,
		vector< double > &	entry
	)

Definition at line 162 of file standaloneTestFastElim.cpp.

Referenced by FastElim::shuffleRows(), and testSorting().

{
    unsigned int num = col.size();
    assert( num == entry.size() );
    // Stupid bubble sort, as we only have up to 5 entries and need to
    // sort both the col and reorder the entries by the same sequence.
    for ( unsigned int i = 0; i < num; ++i ) {
        for ( unsigned int j = 1; j < num; ++j ) {
            if ( col[j] < col[j-1] ) {
                unsigned int temp = col[j];
                col[j] = col[j-1];
                col[j-1] = temp;
                double v = entry[j];
                entry[j] = entry[j-1];
                entry[j-1] = v;
            }
        }
    }
}

Here is the caller graph for this function:

void testSorting

(

)

Definition at line 641 of file standaloneTestFastElim.cpp.

References sortByColumn().

{
    static unsigned int k[] = {20,40,60,80,100,10,30,50,70,90};
    static double d[] = {1,2,3,4,5,6,7,8,9,10};
    vector< unsigned int > col;
    col.insert( col.begin(), k, k+10);
    vector< double > entry;
    entry.insert( entry.begin(), d, d+10);
    sortByColumn( col, entry );
    cout << "testing sorting\n";
    for ( int i = 0; i < col.size(); ++i ) {
        cout << "d[" << i << "]=    " << k[i] <<
           ", col[" << i << "]= " <<    col[i] << ", e=" << entry[i] << endl;
    }
    cout << endl;
}

Here is the call graph for this function:

Variable Documentation

const unsigned int SM_MAX_COLUMNS = 200000

Definition at line 17 of file standaloneTestFastElim.cpp.

Referenced by SparseMatrix< unsigned int >::setSize().

const unsigned int SM_MAX_ROWS = 200000

Template for specialized SparseMatrix. Used both for the Kinetic solver and for handling certain kinds of messages. Speciality is that it can extract entire rows efficiently, for marching through a specified row for a matrix multiplication or for traversing messages.

Requires that type T have an equality operator ==

Definition at line 16 of file standaloneTestFastElim.cpp.

Referenced by SparseMatrix< unsigned int >::setSize().

const unsigned int SM_RESERVE = 8

Definition at line 18 of file standaloneTestFastElim.cpp.

Referenced by SparseMatrix< unsigned int >::SparseMatrix().