home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / SLVBLK.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  72 lines
  1. function x=slvblk(blokmat,b)
  2. % SLVBLK    Solve an almost block-diagonal linear system.
  3. %
  4. %        x=slvblk(blokmat,b)
  5. %
  6. % returns the solution (if any) of the linear system  Ax=b, with the matrix 
  7. %  A  stored in  blokmat  in the spline almost block diagonal form.
  8. %
  9. %  If the system is overdetermined (i.e., has more equations than unknowns),
  10. % the least-squares solution is returned.
  11.  
  12. % C. de Boor / latest change: Jan.14, 1990
  13. % C. de Boor / latest change: Nov.29, 1990 (correct signum problem)
  14. % C. de Boor / latest change: Feb. 1, 1991 (correct signum problem further)
  15. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  16.  
  17.  
  18. % get the basic information
  19. [nb,rows,ncols,last,blocks]=bkbrk(blokmat);
  20.  
  21. ne=sum(rows);nu=sum(last);
  22. if any(cumsum(rows)<cumsum(last))|any(last>ncols),
  23.    error('the coefficient matrix has a nontrivial nullspace')
  24. end
  25.  
  26. [brow,bcol]=size(b);
  27. if(ne~=brow),
  28.    error('matrix and right side are incompatible')
  29. end
  30.  
  31. blocks=[blocks b];
  32. ccols=ncols+bcol;
  33.  
  34. f=1;l=0;
  35. for j=1:nb,
  36.    if (f<=l), % shift the rows still remaining from previous block
  37.       blocks(f:l,:)= ...
  38.          [blocks(f:l,elim+1:ncols) zeros(l+1-f,elim),blocks(f:l,ncols+1:ccols)];
  39.    end
  40.    l=l+rows(j);
  41.    
  42.    elim=last(j);
  43.    % ideally, one would now use
  44.    %   [q,r]=qr(blocks(f:l,1:elim));
  45.    % followed up by
  46.    %   blocks(f:l,:)=q'*blocks(f:l,:);
  47.    %   f=f+elim;
  48.    % but, unfortunately, this generates the possibly very large square matrix q
  49.    % The unhappy alternative is to do the elimination explicitly here, using
  50.    % Householder reflections (and an additional inner loop):
  51.    for k=1:elim,
  52.       a=norm(blocks(f:l,k));
  53.       vv=abs(blocks(f,k))+a; 
  54.       c=vv*a;
  55.       if blocks(f,k)<0, vv=-vv; end
  56.       q=[vv;blocks(f+1:l,k)];
  57.       blocks(f:l,:)=...
  58.        blocks(f:l,:)-((q/c)*ones(1,ccols)).*(ones(l+1-f,1)*(q'*blocks(f:l,:)));
  59.       f=f+1;
  60.    end
  61. end
  62.  
  63. % now we are ready for back-substitution
  64. x=zeros(f-elim-1+ncols,bcol);
  65.  
  66. for j=nb:-1:1,
  67.    elim=last(j);l=f-1;f=f-elim;
  68.    [-x(f-1+[elim+1:ncols],:);eye(bcol)];
  69.    x(f:l,:)=blocks(f:l,1:elim)\(blocks(f:l,elim+1:ccols)*ans);
  70. end
  71. x=x(1:nu,:);
  72.