home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPARFUN.DI$ / SPDIAGS.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.9 KB  |  146 lines
  1. function [res1,res2] = spdiags(arg1,arg2,arg3,arg4)
  2. % SPDIAGS Extract and create sparse band and diagonal matrices.
  3. %
  4. %    SPDIAGS, which generalizes the builtin function "diag", deals with
  5. %    three matrices, in various combinations, as both input and output.
  6. %    A is an m-by-n matrix, usually (but not necessarily) sparse,
  7. %    with its nonzero, or specified, elements located on p diagonals.
  8. %    B is a min(m,n)-by-p matrix, usually (but not necessarily) full,
  9. %    whose columns are the diagonals of A.
  10. %    d is a vector of length p whose integer components specify
  11. %    the diagonals in A.
  12. %
  13. %    Roughly, A, B and d are related by
  14. %        for k = 1:p
  15. %            B(:,k) = diag(A,d(k))
  16. %        end
  17. %
  18. %    Four different operations, distinguished by the number of
  19. %    input arguments, are possible with SPDIAGS:
  20. %
  21. %    Extract all nonzero diagonals:
  22. %        [B,d] = spdiags(A);
  23. %    Extract specified diagonals:
  24. %        B = spdiags(A,d);
  25. %    Replace specified diagonals:
  26. %        A = spdiags(B,d,A);
  27. %    Create a sparse matrix from its diagonals:
  28. %        A = spdiags(B,d,m,n);
  29. %
  30. %    The precise relationship among A, B and d is:
  31. %       if m >= n
  32. %          for k = 1:p
  33. %             for j = max(1,1+d(k)):min(n,m+d(k))
  34. %                B(j,k) = A(j-d(k),j);
  35. %             end
  36. %          end  
  37. %       if m < n
  38. %          for k = 1:p
  39. %             for i = max(1,1-d(k)):min(m,n-d(k))
  40. %                B(i,k) = A(i,i+d(k));
  41. %             end
  42. %          end  
  43. %       end
  44. %    Some elements of B, corresponding to positions "outside" of A,
  45. %    are not defined by these loops.  They are not referenced when
  46. %    B is input and are set to zero when B is output.
  47. %
  48. %    For example, this generates a sparse tridiagonal representation
  49. %    of the classic second difference operator on n points.
  50. %
  51. %        e = ones(n,1);
  52. %        A = spdiags([e -2*e e], -1:1, n, n)
  53. %
  54. %    This turns it into Wilkinson's test matrix (see WILKINSON(n)).
  55. %
  56. %        A = spdiags(abs(-(n-1)/2:(n-1)/2)',0,A)
  57. %
  58. %    Finally, this recovers the 3 diagonals.
  59. %
  60. %        B = spdiags(A)
  61. %
  62. %    The second example is not square.
  63. %
  64. %        A = [ 11    0   13    0
  65. %               0   22    0   24
  66. %               0    0   33    0
  67. %              41    0    0   44
  68. %               0   52    0    0
  69. %               0    0   63    0
  70. %               0    0    0   74]
  71. %
  72. %    has m = 7, n = 4 and p = 3.
  73. %
  74. %    The statement [B,d] = spdiags(A)  produces d = [-3 0 2]'  and
  75. %     
  76. %        B = [ 41   11    0
  77. %              52   22    0
  78. %              63   33   13
  79. %              74   44   24 ]
  80. %
  81. %    Conversely, with the above B and d, the expression
  82. %    spdiags(B,d,7,4) reproduces the original A.
  83. %
  84. %    See also DIAG.
  85.  
  86. %    Rob Schreiber, RIACS, and Cleve Moler, 2/13/91, 6/1/92.
  87. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  88.  
  89.  
  90. if nargin <= 2
  91.    % Extract diagonals
  92.    A = arg1;
  93.    if nargin == 1
  94.       % Find all nonzero diagonals
  95.       [i,j] = find(A);
  96.       d = sort(j-i);
  97.       d = d(find(diff([-inf; d])));
  98.    else
  99.       % Diagonals are specified
  100.       d = arg2;
  101.    end
  102.    [m,n] = size(A);
  103.    p = length(d);
  104.    B = zeros(min(m,n),p);
  105.    for k = 1:p
  106.       if m >= n
  107.          i = max(1,1+d(k)):min(n,m+d(k));
  108.       else
  109.          i = max(1,1-d(k)):min(m,n-d(k));
  110.       end
  111.       B(i,k) = diag(A,d(k));
  112.    end
  113.    res1 = B;
  114.    res2 = d;
  115. end
  116.  
  117. if nargin >= 3
  118.    % Create a sparse matrix from its diagonals
  119.    B = arg1;
  120.    d = arg2;
  121.    p = length(d);
  122.    moda = (nargin == 3);
  123.    if moda
  124.       % Modify a given matrix
  125.       A = arg3;
  126.    else
  127.       % Start from scratch
  128.       A = sparse(arg3,arg4);
  129.    end
  130.    % Process A in compact form
  131.    [i,j,a] = find(A);
  132.    a = [i j a];
  133.    [m,n] = size(A);
  134.    for k = 1:p
  135.       if moda
  136.          % Delete current d(k)-th diagonal
  137.          i = find(a(:,2)-a(:,1) == d(k));
  138.          a(i,:) = [];
  139.       end 
  140.       % Append new d(k)-th diagonal to compact form
  141.       i = (max(1,1-d(k)):min(m,n-d(k)))';
  142.       a = [a; i i+d(k) B(i+(m>=n)*d(k),k)];
  143.    end
  144.    res1 = sparse(a(:,1),a(:,2),full(a(:,3)),m,n);
  145. end
  146.