[Master Index] [Index for PublicToolbox/regutools]

# lanc_b

## (PublicToolbox/regutools/lanc_b.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    [U,B_k,V] = lanc_b(A,p,k,reorth)

### Help Text

    LANC_B Lanczos bidiagonalization.

     B_k = lanc_b(A,p,k,reorth)
     [U,B_k,V] = lanc_b(A,p,k,reorth)

     Performs k steps of the Lanczos bidiagonalization process
     with starting vector p, producing a lower bidiagonal matrix
             [b_11               ]                   [b_21    b_11]
             [b_21 b_22          ]                   [b_32    b_22]
         B = [     b_32 .        ]  stored as  B_k = [ .       .  ] .
             [          . b_kk   ]                   [b_k+1,k b_kk]
             [            b_k+1,k]
     U and V consist of the left and right Lanczos vectors.

     Reorthogonalization is controlled by means of reorth:
        reorth = 0 : no reorthogonalization,
        reorth = 1 : reorthogonalization by means of MGS,
        reorth = 2 : Householder-reorthogonalization.
     No reorthogonalization is assumed if reorth is not specified.

### Cross-Reference Information

This function calls

    app_hh                             C:\BrainStorm_2001\PublicToolbox\regutools\app_hh.m
    gen_hh                             C:\BrainStorm_2001\PublicToolbox\regutools\gen_hh.m

### Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\lanc_b.m

    function [U,B_k,V] = lanc_b(A,p,k,reorth)
    %LANC_B Lanczos bidiagonalization.
    %
    % B_k = lanc_b(A,p,k,reorth)
    % [U,B_k,V] = lanc_b(A,p,k,reorth)
    %
    % Performs k steps of the Lanczos bidiagonalization process
    % with starting vector p, producing a lower bidiagonal matrix
    %         [b_11               ]                   [b_21    b_11]
    %         [b_21 b_22          ]                   [b_32    b_22]
    %     B = [     b_32 .        ]  stored as  B_k = [ .       .  ] .
    %         [          . b_kk   ]                   [b_k+1,k b_kk]
    %         [            b_k+1,k]
    % U and V consist of the left and right Lanczos vectors.
    %
    % Reorthogonalization is controlled by means of reorth:
    %    reorth = 0 : no reorthogonalization,
    %    reorth = 1 : reorthogonalization by means of MGS,
    %    reorth = 2 : Householder-reorthogonalization.
    % No reorthogonalization is assumed if reorth is not specified.

    % Reference: G. H. Golub & C. F. Van Loan, "Matrix Computations",
    % 2. Ed., Johns Hopkins, 1989.  Section 9.3.4.

    % Per Christian Hansen, UNI-C, 05/25/93.

    % Initialization.
    if (k<1), error('Number of steps k must be positive'), end
    if (nargin < 4), reorth = 0; end
    if (reorth < 0 | reorth > 2), error('Illegal reorth'), end
    if (nargout==2), error('Not enough output arguments'), end
    [m,n] = size(A);
    if (nargout>1 | reorth==1)
      U = zeros(m,k); V = zeros(n,k); UV = 1;
    else
      UV = 0;
    end
    if (reorth==2)
      if (k>=n), error('No. of iterations must satisfy k < n'), end
      HHU = zeros(m,k); HHV = zeros(n,k);
      HHalpha = zeros(1,k); HHbeta = HHalpha;
    end

    % Prepare for Lanczos iteration.
    v = zeros(n,1);
    beta = norm(p);
    if (beta==0), error('Starting vector must be nonzero'), end
    if (reorth==2)
      [beta,HHbeta(1),HHU(:,1)] = gen_hh(p);
    end
    u = p/beta;
    if (UV), U(:,1) = u; end

    % Perform Lanczos bidiagonalization with/without reorthogonalization.
    for i=1:k

      r = A'*u - beta*v;
      if (reorth==0)
        alpha = norm(r); v = r/alpha;
      elseif (reorth==1)
        for j=1:i-1, r = r - (V(:,j)'*r)*V(:,j); end
        alpha = norm(r); v = r/alpha;
      else
        for j=1:i-1
          r(j:n) = app_hh(r(j:n),HHalpha(j),HHV(j:n,j));
        end
        [alpha,HHalpha(i),HHV(i:n,i)] = gen_hh(r(i:n));
        v = zeros(n,1); v(i) = 1;
        for j=i:-1:1
          v(j:n) = app_hh(v(j:n),HHalpha(j),HHV(j:n,j));
        end
      end
      B_k(i,2) = alpha; if (UV), V(:,i) = v; end

      p = A*v - alpha*u;
      if (reorth==0)
        beta = norm(p); u = p/beta;
      elseif (reorth==1)
        for j=1:i, p = p - (U(:,j)'*p)*U(:,j); end
        beta = norm(p); u = p/beta;
      else
        for j=1:i
          p(j:m) = app_hh(p(j:m),HHbeta(j),HHU(j:m,j));
        end
        [beta,HHbeta(i+1),HHU(i+1:m,i+1)] = gen_hh(p(i+1:m));
        u = zeros(m,1); u(i+1) = 1;
        for j=i+1:-1:1
          u(j:m) = app_hh(u(j:m),HHbeta(j),HHU(j:m,j));
        end
      end
      B_k(i,1) = beta; if (UV), U(:,i+1) = u; end

    end

    if (nargout==1), U = B_k; end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON