|
Boost Users : |
From: nisha kannookadan (nishak44_at_[hidden])
Date: 2007-11-30 12:05:55
Hi guys
Thanks for all. I follow your advises, posted the code in the ublas, gonna read a C++ book
and Im positive, its gonna be fine.
I changed now the code to something runnable and have it also in matlab,
u will see, its way slower in C++.
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
But I really gotta say, this mailing list is great. I really appreciate all of your replies.
Thanks,
Nisha K
_________________________________________________________________
News, entertainment and everything you care about at Live.com. Get it now!
http://www.live.com/getstarted.aspx
Boost-users list run by williamkempf at hotmail.com, kalb at libertysoft.com, bjorn.karlsson at readsoft.com, gregod at cs.rpi.edu, wekempf at cox.net