merge with beta

ABfgPreEliminate.m

@@ -1,96 +1,171 @@
-function [x,y,tolA,tolB]=ABfgPreEliminate(CtrlVar,A,B,f,g)
-
-
-
-    [nA,mA]=size(A);
-    [nB,mB]=size(B);
-    [nf,mf]=size(f);
-
-    if isempty(B) && isempty(g) && ~isempty(A) && ~isempty(f) && mA==nf
-
-        % Possibly not needed, but check if this is not just a very simple case of B=g=[]
-
-        x=A\f;
-        y=NaN;
-
-        if nargout>2
-            tolB=NaN;
-            tolA=norm(A*x-f)/norm(f);
+function [x,y,dAtilde,tolA,tolB]=ABfgPreEliminate(CtrlVar,A,B,f,g,dAtilde)
+
+
+narginchk(5,6)
+
+if nargin< 6
+    dAtilde=[];
+end
+
+[nA,mA]=size(A);
+[nB,mB]=size(B);
+[nf,mf]=size(f);
+
+
+if isempty(B) && isempty(g) && ~isempty(A) && ~isempty(f) && mA==nf
+
+    % Possibly not needed, but check if this is not just a very simple case of B=g=[]
+
+    x=A\f;
+    y=NaN;
+
+    if nargout>3
+        tolB=NaN;
+        tolA=norm(A*x-f)/norm(f);
+    end
+
+else
+
+    BBT=B*B';
+
+    if isdiag(BBT)  % the method assumes that B B' is diagonal
+
+
+        % It also assumes that B B' is a unity matrix, but if not then simple scaling can be used
+        % to ensure that this is the case.
+        % To make this a bit more general, I here check if B B' is indeed unity, and
+        % if not I do the requried scaling.
+        tolerance=eps*1000;
+        isBBTunity=all(abs(diag(BBT) - 1) < tolerance) ;
+
+        if ~isBBTunity
+            [B,g,~,Scale]=ScaleL(CtrlVar,B,g) ;
+        else
+            Scale=1;
         end
+
+        % For numerical reasons a further simple scaling of A is done to bring the
+        % sizes of the elements of A in line with those of B.
+        factor=1./(full(mean(abs(diag(A)))));
+
+        if ~isfinite(factor)  % just in case all elements along the diagonal happen to be equal to zero
+            factor=1;
+        end
+
+        BtB=B'*B;
+
+        A=factor*A ; f=factor*f ;  % this leaves x unaffected but y is scaled
+
+
+
+        Q=speye(nA,nA)-BtB ;
+        Atilde=Q*A+ BtB ;
+        btilde=(Q*f+B'*g) ;
+
+
+        if CtrlVar.Parallel.isTest
+
+            % seq:
+
+            % distribute
+            tSeq=tic ;
+            dAtildeSeq=decomposition(Atilde);  % this might not be needed if this has been done already, but for speed comparison this is done here each time
+            xSeq=dAtildeSeq\btilde;
+            tSeq=toc(tSeq) ; 
+
+
+            % distribute
+            tDistributed=tic ;
+            AtildeDist=distributed(Atilde);
+            btildeDist=distributed(btilde);
+            dAtildeDist=decomposition(AtildeDist);  % this might not be needed if this has been done already, but for speed comparison this is done here each time
+            xDist=dAtildeDist\btildeDist;
+            tDistributed=toc(tDistributed) ; 
 
-    else
-
-        BBT=B*B';
-
-        if isdiag(BBT)  % the method assumes that B B' is diagonal
-
-
-            % It also assumes that B B' is a unity matrix, but if not then simple scaling can be used
-            % to ensure that this is the case.
-            % To make this a bit more general, I here check if B B' is indeed unity, and
-            % if not I do the requried scaling.
-            tolerance=eps*1000;
-            isBBTunity=all(abs(diag(BBT) - 1) < tolerance) ;
-
-            if ~isBBTunity
-                [B,g,~,Scale]=ScaleL(CtrlVar,B,g) ;
-            else
-                Scale=1;
+            [nAtilde,mAtilde]=size(Atilde);
+            fprintf('\n ----------------------------- Info on distributed solve performance : \n')
+            fprintf('%i x %i :  tSeq=%f \t tDistributed=%f \t tSeq/rDistributed=%f \n',nAtilde,mAtilde,tSeq,tDistributed,tSeq/tDistributed) ;
+            fprintf(' norm(xSeq-xDist)/norm(xSeq)=%g    \n',full(norm(xSeq-xDist)/norm(xSeq)))
+            fprintf(' ----------------------------- \n')
+
+
+        end
+
+
+        % https://uk.mathworks.com/help/parallel-computing/benchmarking-a-b.html
+        if CtrlVar.Parallel.Distribute
+            if ~isdistributed(Atilde)
+                Atilde=distributed(Atilde);
             end
-
-            % For numerical reasons a further simple scaling of A is done to bring the
-            % sizes of the elements of A in line with those of B.
-            factor=1./(full(mean(abs(diag(A)))));
-
-            if ~isfinite(factor)  % just in case all elements along the diagonal happen to be equal to zero
-                factor=1;
+            if ~isdistributed(btilde)
+                btilde=distributed(btilde);
             end
-
-            BtB=B'*B;
-
-            A=factor*A ; f=factor*f ;  % this leaves x unaffected but y is scaled
-
-
-
-            Q=speye(nA,nA)-BtB ; 
-            Atilde=Q*A+ BtB ;
-            btilde=(Q*f+B'*g) ;
-
-
-            %dAtilde=factorize(Atilde);  % for some reason, this does make things slower
-
-            x=Atilde\btilde;
-            y=B*(f-A*x);
-
-            % Now the solution of the scaled system has been found.
-
-            y=y/factor;
-
-            if nargout>2
-                % check if within tolerances
-
-                A=A/factor ;
-                f=f/factor ;
-                tolA=norm(A*x+B'*y-f)/norm(f);
-                tolB=norm(B*x-g);
-
-                if tolA>1e-6 || tolB>1e-6
-
-                    fprintf('ABfgPreEliminate: Error seems too large or \t \t \t %g \t %g \n ',norm(A*x+B'*y-f)/norm(f),norm(B*x-g))
-
-                end
+        end
+
+        %  decomposition is about the same, and as expected this only speeds things up if several solves with the same matrix
+        %  are needed.
+        %
+        tDecomposition=tic;
+        if isempty(dAtilde)
+            dAtilde=decomposition(Atilde);
+        end
+        tDecomposition=toc(tDecomposition) ;
+
+        tSolve=tic;
+        x=dAtilde\btilde;
+        tSolve=toc(tSolve);
+
+        % tSolve=tic;
+        % x=Atilde\btilde;
+        % tSolve=toc(tSolve);
+
+        % [tDecomposition tSolve]
+
+        if isdistributed(x)
+            x=gather(x) ;
+        end
+        % toc
+
+
+
+
+
+
+
+
+
+
+        y=B*(f-A*x);
+
+        % Now the solution of the scaled system has been found.
+
+        y=y/factor;
+
+        if nargout>3
+            % check if within tolerances
+
+            A=A/factor ;
+            f=f/factor ;
+            tolA=norm(A*x+B'*y-f)/norm(f);
+            tolB=norm(B*x-g);
+
+            if tolA>1e-6 || tolB>1e-6
+
+                fprintf('ABfgPreEliminate: Error seems too large or \t \t \t %g \t %g \n ',norm(A*x+B'*y-f)/norm(f),norm(B*x-g))
+
             end
-
-            y=Scale*y; % and now scale y in case B and g were scaled above.
-
-        else
-
-
-            error('ABfgPreEliminate:B','B*B^T not diagonal')
         end
-
+
+        y=Scale*y; % and now scale y in case B and g were scaled above.
+
+    else
+
+
+        error('ABfgPreEliminate:B','B*B^T not diagonal')
     end
-
+
+end
+
 end