[ublas] C++ Performance

From: nisha kannookadan (nishak44_at_[hidden])
Date: 2007-11-30 12:03:30

Hey everybody

Finally I got to the ublas mailing list and I hope u guys can help me.
Im quite new with boost and everything else than an expert.

Im programming in C++ and have to translate a program from matlab
to C++. I did it, but my program is way to slow.

I got already suggestions from the boost users and developpers:
- to do pass by refernce
- define DNDEBUG
- to compile with full optimization
- to use bindings (but I dont think, I need it, since its not high level matrix stuff)
- and that I should go to post the question here :).

I did some of the changes already and it got little better, but its still bad.

The program is quite huge, so I took a part out of it, which is now even runnable
for you.
I post here the maltab and c++ code, so u maybe can help me.

So first C++:

#include 
#include "MVHF.h"
#include 
#include 
#include 
#include 
#include 

#define DNDEBUG
#include 

using namespace std;

void ttrans(matrix& At, int level) 
{ 
        matrix  cfe1, cfe2, cfo, cfe, c, d;
        int N,s2;
        
        N = (At.size1()+1)/2;
        s2 = At.size2();
        zero_matrix zer(N,s2);
        
        for (int ii = 1; ii <= level; ii++)
        {
                
                cfo = subslice(At, 0,2,N, 0,1,s2);
                cfe = subslice(At, 1,2,N-1, 0,1,s2);
                
                c = cfe + (subrange(cfo, 0,N-1, 0,s2)+subrange(cfo, 1,N, 0,s2))*0.5;
                
                zer.resize(N,s2,true);
                cfe1 = zer;
                cfe2 = zer;
                
                (subrange(cfe1, 0,N-1, 0,s2)).assign(cfe);
                (subrange(cfe2, 1,N, 0,s2)).assign(cfe);
                d = cfo-(cfe1+cfe2)*0.5;
                
            (subrange(At, 0,N-1, 0,At.size2())).assign(c);
            (subrange(At, N-1,2*N-1, 0,At.size2())).assign(d);
            
            N = N/2;
        }
        
}

void init(int dom,int llev)
{	
                        
        int in_dofs_LA= 511;
        matrix Am_LA, As_LA;
        double mesh_size = ((double)(dom)/(in_dofs_LA+1));
                                
        scalar_matrix zer(in_dofs_LA, in_dofs_LA,0.0);
        Am_LA.resize(in_dofs_LA, in_dofs_LA, false);
        As_LA.resize(in_dofs_LA, in_dofs_LA, false);
        Am_LA.assign(zer);
        As_LA.assign(zer);
                                
        Am_LA(0,0) = 4.0;
        Am_LA(1,0) = 1.0;
                                
        As_LA(0,0) = 2.0;
        As_LA(1,0) = -1.0;
                                
        for (int ii = 1; ii < (in_dofs_LA-1); ii++){
                Am_LA(ii-1,ii) = 1.0;
                Am_LA(ii,ii) = 4.0;
                Am_LA(ii+1,ii) = 1.0;
                                        
                As_LA(ii-1,ii) = -1.0;
                As_LA(ii,ii) = 2.0;
                As_LA(ii+1,ii) = -1.0;	
        }
                                
        Am_LA(in_dofs_LA-2,in_dofs_LA-1) = 1.0;
        Am_LA(in_dofs_LA-1,in_dofs_LA-1) = 4.0;
                                
        As_LA(in_dofs_LA-2,in_dofs_LA-1) = -1.0;	
        As_LA(in_dofs_LA-1,in_dofs_LA-1) = 2.0;
                                
        Am_LA = Am_LA*(mesh_size/6.0);
        As_LA = As_LA*(1.0/mesh_size);
        
        ttrans(Am_LA,llev);
        Am_LA = trans(Am_LA);
        ttrans(Am_LA,llev);
        Am_LA = trans(Am_LA);
                        
        ttrans(As_LA,llev);
        As_LA = trans(As_LA);
        ttrans(As_LA,llev);
        As_LA = trans(As_LA);
        
}

int main() {
        clock_t start, end;
        start = clock();

        init(1,8);
        end = clock();
        
        std::cout << " " << std::endl;
        std::cout << "Elapsed time is " << (double)(end-start)/CLOCKS_PER_SEC << "." << std::endl;
}

Here the matlab version of it:

function test()
tic
R = 1;   
L = 8;
n = 2^(L+1)-1;                            
h = R/(n+1);                             

e = ones(n,1);
Am = h/6*spdiags([e, 4*e, e], -1:1, n, n);

% compute stiffness matrix
As = 1/h*spdiags([-e, 2*e, -e], -1:1, n, n);

% transform into wavelets
Am = ttrans(ttrans(Am,L)',L)';
As = ttrans(ttrans(As,L)',L)';
toc

return

function cd = ttrans(cl,L)
N = size(cl,1)+1;
N = N/2;

m = size(cl,2);

z = zeros(1,m);
cd = cl;
for i=1:L
    cfo = cd(1:2:2*N-1,:);
    cfe = cd(2:2:2*N-2,:);
    
    % wavelets
    c = cfe + (cfo(1:end-1,:)+cfo(2:end,:))/2;
    d = cfo - ([cfe;z]+[z;cfe])/2;
     
    cd(1:2*N-1,:) = [c;d];
    N = N/2;
end

return

Id be so thankful for any help. 
Bye
Nisha K
_________________________________________________________________
Invite your mail contacts to join your friends list with Windows Live Spaces. It's easy!
http://spaces.live.com/spacesapi.aspx?wx_action=create&wx_url=/friends.aspx&mkt=en-us