spm_kron.m 1.2 KB

12345678910111213141516171819202122232425262728293031323334353637383940
  1. function K = spm_kron(A,B)
  2. % as for kron, but forces sparse outputs
  3. %KRON Kronecker tensor product.
  4. % KRON(X,Y) is the Kronecker tensor product of X and Y.
  5. % The result is a large matrix formed by taking all possible
  6. % products between the elements of X and those of Y. For
  7. % example, if X is 2 by 3, then KRON(X,Y) is
  8. %
  9. % [ X(1,1)*Y X(1,2)*Y X(1,3)*Y
  10. % X(2,1)*Y X(2,2)*Y X(2,3)*Y ]
  11. %
  12. % Class support for inputs X,Y:
  13. % float: double, single
  14. % Previous versions by Paul Fackler, North Carolina State,
  15. % and Jordan Rosenthal, Georgia Tech.
  16. % Copyright 1984-2004 The MathWorks, Inc.
  17. % $Revision: 2863 $ $Date: 2004/06/25 18:52:18 $
  18. %__________________________________________________________________________
  19. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  20. % Karl Friston
  21. % $Id: spm_kron.m 2863 2009-03-11 20:25:33Z guillaume $
  22. [ma,na] = size(A);
  23. [mb,nb] = size(B);
  24. [ia,ja,sa] = find(A);
  25. [ib,jb,sb] = find(B);
  26. ia = ia(:); ja = ja(:); sa = sa(:);
  27. ib = ib(:); jb = jb(:); sb = sb(:);
  28. ka = ones(size(sa));
  29. kb = ones(size(sb));
  30. t = mb*(ia-1)';
  31. ik = t(kb,:)+ib(:,ka);
  32. t = nb*(ja-1)';
  33. jk = t(kb,:)+jb(:,ka);
  34. K = sparse(ik,jk,sb*sa.',ma*mb,na*nb);