opSubsAsgn.m: Fixed bug in scalar assignment

opBlockDiag.m: Fixed bug with weights double.m: Use minimal number of multipliations git-svn-id: svn://tuatara.cs.ubc.ca/spot/trunk@18 602325e8-326d-0410-a6db-8800ab266cdc
mpf · Jun 30, 2009 · 8a771cc · 8a771cc
1 parent 7edcd3a
commit 8a771cc
Show file tree

Hide file tree

Showing 6 changed files with 38 additions and 25 deletions.
diff --git a/@opSpot/diag.m b/@opSpot/diag.m
@@ -0,0 +1,15 @@
+function d = diag(A)
+
+%   Copyright 2009, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: double.m 13 2009-06-28 02:56:46Z mpf $
+
+[m,n] = size(A);
+k = min(m,n);
+
+d = zeros(k,1);
+for i=1:k
+   v = zeros(n,1); v(i) = 1;
+   w = A*v;
+   d(i) = w(i);
+end
diff --git a/@opSpot/double.m b/@opSpot/double.m
@@ -8,4 +8,8 @@
 %   http://www.cs.ubc.ca/labs/scl/sparco
 %   $Id$
 
-M = A*speye(size(A,2));
+if size(A,1) < size(A,2)
+   M = (A'*speye(size(A,1)))';
+else
+   M = A*speye(size(A,2));
+end
diff --git a/@opSpot/subsasgn.m b/@opSpot/subsasgn.m
@@ -114,7 +114,7 @@
       % ---------------------------------------------------------------
       % Check if size of indices and b match
       % ---------------------------------------------------------------
-      if ~(size(b,1) == sizeIndex(1)) || ~(size(b,2) == sizeIndex(2))
+      if ~isscalar(b) && (size(b,1) ~= sizeIndex(1) || size(b,2) ~= sizeIndex(2))
          error('Subscripted assignment dimension mismatch.');
       end
 

diff --git a/TODO b/TODO
@@ -1,18 +1,13 @@
 TODO
-[] opBlockDiag.m; Is the parameter input okay?
-                  Avoid replicating operator explicitly?
+[] opBlockDiag.m; Avoid replicating operator explicitly?
                   How to display repeated blocks: opBlockDiag(ones(10,1),opDCT(4))?
-   DECISION: Fix weights, add to multiplication; make sure help
-   comment mentions weight vector in combination with a single
+   make sure help comment mentions weight vector in combination with a single
    operator as well.
-[] @opSpot/diag.m; Should this be similar to blkdiag? NO, just like
-   matlab's reverse diag.
+[] @opSpot/diag.m; add documentation
 
 [] opInverse.m: Should we allow operators such as opInverse on
    non-linear operators? NO!!! Check linear flag. LIKEWISE for PCG, LSQR, ...
 
-
-[] @opSpot/double.m; Check dimensions and minimize number of matrix-vector products
 [] Clean up multiplication implementations, see e.g., opCurvelet
 [] opWavelet.m: Help comment
 [] opToeplitz.m, opToepGauss.m, opToepSign.m: Help comment

diff --git a/opBlockDiag.m b/opBlockDiag.m
@@ -125,7 +125,7 @@
 
           % Construct function handle
           if overlap == 0
-             fun = @(x,mode) opBlockDiag_intrnl(m,n,opList,x,mode);
+             fun = @(x,mode) opBlockDiag_intrnl(m,n,opList,weights,x,mode);
           elseif overlap < 0
              % Overlap in columns
              m = 0; n = 0; colOffset = 0; column = 0;
@@ -139,7 +139,7 @@
                 column    = column + nOp + overlap; % Overlap is negative
              end
              n = n - colOffset;
-             fun = @(x,mode) opBlockDiagCol_intrnl(m,n,-colOffset,-overlap, opList,x,mode);
+             fun = @(x,mode) opBlockDiagCol_intrnl(m,n,-colOffset,-overlap,opList,weights,x,mode);
           else
              % Overlap in rows
              m = 0; n = 0; rowOffset = 0; row = 0;
@@ -153,7 +153,7 @@
                 row       = row + mOp - overlap; % Overlap is positive
              end
              m = m - rowOffset;
-             fun = @(x,mode) opBlockDiagRow_intrnl(m,n,-rowOffset,overlap,opList,x,mode);
+             fun = @(x,mode) opBlockDiagRow_intrnl(m,n,-rowOffset,overlap,opList,weights,x,mode);
           end
 
           % Construct operator
@@ -203,7 +203,7 @@
 %=======================================================================
 
 
-function y = opBlockDiag_intrnl(m,n,opList,x,mode)
+function y = opBlockDiag_intrnl(m,n,opList,weights,x,mode)
 
 kx = 0; ky = 0;
 
@@ -212,7 +212,7 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+mOp) = op * x(kx+1:kx+nOp);
+      y(ky+1:ky+mOp) = weights(i) * op * x(kx+1:kx+nOp);
       kx = kx + nOp;
       ky = ky + mOp;
    end;
@@ -221,7 +221,7 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+nOp) = op' * x(kx+1:kx+mOp);
+      y(ky+1:ky+nOp) = conj(weights(i)) * op' * x(kx+1:kx+mOp);
       kx = kx + mOp;
       ky = ky + nOp;
    end
@@ -232,7 +232,7 @@
 %=======================================================================
 
 
-function y = opBlockDiagCol_intrnl(m,n,offset,overlap,opList,x,mode)
+function y = opBlockDiagCol_intrnl(m,n,offset,overlap,opList,weights,x,mode)
 
 kx = 0; ky = 0;
 
@@ -241,7 +241,7 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+mOp) = op * x(kx+1:kx+nOp);
+      y(ky+1:ky+mOp) = weights(i) * (op * x(kx+1:kx+nOp));
       kx = kx + nOp - overlap;
       ky = ky + mOp;
    end
@@ -250,7 +250,7 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+nOp) = y(ky+1:ky+nOp) + op' * x(kx+1:kx+mOp);
+      y(ky+1:ky+nOp) = y(ky+1:ky+nOp) + conj(weights(i))*(op' * x(kx+1:kx+mOp));
       kx = kx + mOp;
       ky = ky + nOp - overlap;
    end   
@@ -260,7 +260,7 @@
 %=======================================================================
 
 
-function y = opBlockDiagRow_intrnl(m,n,offset,overlap,opList,x,mode)
+function y = opBlockDiagRow_intrnl(m,n,offset,overlap,opList,weights,x,mode)
 
 kx = 0; ky = 0;
 
@@ -269,7 +269,7 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+mOp) = y(ky+1:ky+mOp) + op * x(kx+1:kx+nOp);
+      y(ky+1:ky+mOp) = y(ky+1:ky+mOp) + weights(i) * (op * x(kx+1:kx+nOp));
       kx = kx + nOp;
       ky = ky + mOp - overlap;
    end
@@ -278,9 +278,9 @@
    for i=1:length(opList)
       op = opList{i};
       [mOp,nOp] = size(op);
-      y(ky+1:ky+nOp) = op' * x(kx+1:kx+mOp);
+      y(ky+1:ky+nOp) = conj(weights(i)) * (op' * x(kx+1:kx+mOp));
       kx = kx + mOp - overlap;
       ky = ky + nOp;
    end   
 end
-end
+end
diff --git a/opSubsAsgn.m b/opSubsAsgn.m
@@ -109,12 +109,11 @@
 
           % Check if index and operator dimensions match
           if isscalar(B)
-             B = opOnes(sizeIndex(1),sizeIndex(2))*double(B);
+             B = opFoG(B,opOnes(sizeIndex(1),sizeIndex(2)));
           elseif (sizeIndex(1) ~= size(B,1)) || (sizeIndex(2) ~= size(B,2))
              error('Subscripted assignment dimension mismatch.');
           end
 
-
           % Determine final operator size
           m = max(maxIndex(1),size(A,1));
           n = max(maxIndex(2),size(A,2));