function v = svec(A,blk,sp_blk)
% SVEC   creates a vector v of length sum_i  blk(i)*(blk(i)+1)/2,
%        (sum over # blocks) from a symmetric nxn matrix A
%
% v = svec(A,blk,sp_blk)
%
%   Here n = sum_i blk(i); svec stacks the columns of the lower
%   triangular part of A on top of one another, after multyplying
%   the off diagonal entries by sqrt(2).
%   The factor sqrt(2) is to ensure that svec is an isometry, i.e.
%                  <M,N> = vec(M)^T vec(N)
%   for any pair of symmetric matrices M, N.
%
%   Note: svec only uses the lower triangular part of the matrix A!
%
% input variables:
%     - A         a block diagonal symmetric matrix
%     - blk       block structure vector
%     - sp_blk    flag: if sp_blk = 1, svec returns a sparse vector and a full vector otherwise
%
% output variables:
%     - v         the vector representation of the block diagonal symmetric matrix

% ??????? Version 0.92 BETA
% Copyright (c) 1997 by
% J.-P. Haeberly, M. Nayakkankuppam, M.L. Overton
% Last modified: 9/3/97
%
 sparseblk = 0;
 if nargin == 3,
    if sp_blk == 1,
       sparseblk = 1;
    end;
 end;
 d = diag(A);          % save the old diagonal
 A = A - diag(d);      % A's diagonal elts are now zero
 A = sqrt(2)*A;
 A = A + diag(d);      % restore diagonal
%
 start = 1; fin = 0;
 vstart = 1; vfin = 0;
 for bidx = 1:length(blk),
    bsize = blk(bidx);
    b2size = bsize*(bsize+1)/2;
    fin = fin + bsize;
    vfin = vfin + b2size;
    block = A(start:fin,start:fin);
    J = (1:bsize)';
    J = J(:,ones(1,bsize));
    I = J';
    if bidx == 1,
       v = block(I<=J);
    else
        v(vstart:vfin,1) = block(I<=J);
    end;
    start = start + bsize;
    vstart = vstart + b2size;
 end;   % bidx loop
%
 if sparseblk,
    v = sparse(v);
 end;
%
% END function
