Re: [ublas] Solving linear equations in a loop (Optimization issues)

From: Sourabh (sourabh_at_[hidden])
Date: 2007-03-03 07:27:56

On Fri, 2 Mar 2007, Gunter Winkler wrote:

> On Friday 02 March 2007 15:12, Karl Meerbergen wrote:
> > It should be possible to do this with the triangular solves from ublas
> > but I am unsure whether they efficiently handle sparse matrices. Gunter,
> > what about the new routines?
> 
> I have specializations for compressed_matrix. They are hopefully as fast as 
> possible ;-)
> 
>   template < bool unit, class T, class Z, class D, size_t IB, class IA, 
>              class TA, class VEC >
>   bool
>   inplace_solve_lower(
>     const compressed_matrix<T, basic_row_major<Z,D>, IB, IA, TA>& L,
>     vector_expression<VEC>& ve)
>   {
> 	...
>   }
> 

The program (see below), after some modifications is still taking around 
8 minutes for the arguments 1 200000 100. 

Even if I remove solve, and addition methods, the statement
 lembda1.assign (ident-lembda) is causing the program to take 1 minute.

I am not sure where I can optimize it better ?
 

int main (int argc, char* argv[]) 
{
    unsigned int samples    = atoi (argv[1]);
    unsigned int numpaths   = atoi (argv[2]);
    unsigned short gatesperpath = atoi (argv[3]);
    identity_matrix<double> ident (gatesperpath);
    vector<double> slopes;
    vector <double> delays (gatesperpath);

    for (unsigned int i = 0; i < samples; ++i) {
        vector <double> omegas(gatesperpath);
        for (unsigned int j = 0; j < numpaths ; ++ j) {
            compressed_matrix<double> 	lembda (gatesperpath, gatesperpath);
            compressed_matrix<double> 	lembda1 (gatesperpath,gatesperpath);
            lembda1.assign (ident-lembda);
            slopes = solve (lembda1, omegas, lower_tag ());
            noalias( delays ) = omegas;
            axpy_prod( lembda1, slopes, delays, false );
        }
        
    }
    
}