Home > source > trdsdsmat.m

# trdsdsmat

## PURPOSE

compute trace(A*S1*B*S2) where A,B are dense symmetric matrices, S1,S2 are sparse

## SYNOPSIS

function [tr]=trdsdsmat(A,S1,B,S2)

## DESCRIPTION

     compute trace(A*S1*B*S2) where A,B are dense symmetric matrices, S1,S2 are sparse
     typically, it would be: trdsdsmat(pZUZ, Akdix, invZ, Akdjx)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 % compute trace(A*S1*B*S2) where A,B are dense symmetric matrices, S1,S2 are sparse
    0002 % typically, it would be: trdsdsmat(pZUZ, Akdix, invZ, Akdjx)
    0003 function [tr]=trdsdsmat(A,S1,B,S2)
    0004
    0005   %tic;
    0006   %S1BS2=S1*B*S2;
    0007   %tr = trace(A*S1BS2);
    0008   %toc
    0009
    0010   %tic;
    0011   %S1BS2=sparse(S1*(full(B)*S2)); % or sparse(full(B)*S2)...?
    0012   %tr = trace(A*S1BS2);
    0013   %toc
    0014
    0015   % USE ME (if not mex)!   % -->7.5s (sparse), ?? (dense)
    0016   %tic;
    0017   S1BS2=S1*(B*S2); % or sparse(full(B)*S2)...?
    0018   tr = A(:)'*S1BS2(:);
    0019   %toc
    0020
    0021   % or use mex
    0022   %tr=mextrdsdsmat(A,S1,full(B),S2);  % very slow... --> 15s
    0023   %%tr=mextrdsdsmat(A,S1,B,S2);  % with full one above --> ~6s
    0024
    0025   %S1BS2=S1*sparse((full(B)*S2)); % or sparse(full(B)*S2)...?
    0026   %tr = A(:)'*S1BS2(:);
    0027
    0028   % rather list the elements...?
    0029   % zkusit tr(A,B) = sum A_ij * B_ij = svec(A)'svec(B)
    0030   % careful about svec on symmetric matrices (sqrt(2)...)
    0031
    0032