Re: R: R: help on Octave

[c.]

On 29 Jun 2012, at 12:33, Pierpaolo wrote:

> Hi, 
> 
> I created the sparse matrix and sparse array by
> 
> s=speye(3*dd);  
> d=speye(3*dd,1);
> 
> This because I don't know how to create an empty sparse matrix:
> w=s-s;   
> jj=s-s;
> fe=d-d;

to allocate an empty matrix you can do 

s = sparse (3*dd, 3*dd);

I see no point in allocating an empty sparse vector.

> Then I fill the matrix by several for loop like this:
> 
> function jac=fja(x)
> jac=2/(1-x)^3;
> endfunction
> 
> for j=1:(dimytot-4) 
>       for i=1:(dimxtot-4) 
>       
> jj((j+1)*dimxtot+(i+2),(j+1)*dimxtot+(i+2))=feval("fja",vp((j+1)*dimxtot+(i+
> 2)));
>       endfor
> endfor
> 
> 
> [L,U]=lu((w-jj));
> 
> err=U\(L\(fe-w*vp));

It would be _MUCH_ more efficient if you would first create your amtrix in AIJ 
format and then convert to a sparse matrix, i.e.

idx = 1;
for j=1:(dimytot-4) 
  for i=1:(dimxtot-4) 
   irow(idx) = (j+1)*dimxtot+(i+2);
   jcol(idx) = (j+1)*dimxtot+(i+2); 
   jj(idx++) = feval("fja",vp((j+1)*dimxtot+(i+2)));
  endfor
endfor

jj = sparse (irow, jcol, jj, 3*dd, 3*dd);

BTW, what's the reason for 

feval("fja",vp((j+1)*dimxtot+(i+2)))

why not just 

fja (vp((j+1)*dimxtot+(i+2)))

?


HTH,
c.