%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  MAP05 - Metodos Matriciais Computacionais - PPGMAp UFRGS
%%  Mat01050 - Algebra Matricial Computacional - DMPA UFRGS
%%
Tecnicamente, esse nao eh um script, pois nao pode ser executado.
Entenda como um conjunto de instrucoes que facilmente podem
ser carregados com o mouse.
Tarefa: implementar interface entre Scilab e rotina dpbtrf (fatoracao
de Cholesky de matriz banda) da biblioteca Lapack

// execute no Scilab
link("/usr/lib/liblapack.so","dpbtrf");
//
function L=lapack_chol_band(A,p)
  // fatorizacao de Cholesky: A deve ser simetrica positiva definida
  // p deve ser o numero de sub-diagonais
  // somente a parte inferior da matriz A eh considerada
  m=size(A,1); n=size(A,2);
  if (m~=n) then error('matriz deve ser quadrada'); end;
  uplo='L'; ab=zeros(p+1,n); ldab=p+1; //
  for j=1:n
    ab(1:min(j+p,n)-j+1,j)= A(j:min(j+p,n),j);
  end
  vl = call('dpbtrf',uplo,1,"c",n,2,"i",p,3,"i",ab,4,"d",ldab,5,"i","out",[ldab,m],4,"d");
  L=0; for j=1:n
      L( j:j+p,j) = vl( :,j );
  end
  L = L(1:n,:);  // remove linhas excedentes na parte inferior
endfunction

//
n=5; A=diag(4*ones(n,1));
for i=1:p
  v=2*rand(n-i,1)-ones(n-i,1);A=A + diag(v,i)+diag(v,-i);
end
//
L = lapack_chol_band(A)
norm(A-L*L')
%
% JBC 21Nov2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%