FileMatrix.m

Go to the documentation of this file.

namespace data{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class FileMatrix
  :public data.FileData,
   public data.ABlockedData {
    public: /* ( Constant ) */

            static const .integer BLOCK_SIZE = 512;
    public:

        InPlaceTranspose = false;
    public:

        integer bCols;
        integer nBlocks;
        integer n;
        integer m;
        MinValue = Inf;


        MaxValue = -Inf;

        
    private:

        idx;
        created;
        blocksize;
    private:

        cachedBlock = "[]";
        cachedNr = 0;

        cacheDirty = false;

        
    public:

        FileMatrix(matrix<double> var,varargin) {

            if isempty(var)
                error(" Cannot create a FileMatrix with empty first argument. Must either be a row number or a matrix. ");
            end

            
            matrixin = false;
            if numel(var) > 1
                A = var;
                [n, m] = size(A); /* #ok<*PROP> */

                matrixin = true;
/*              elseif isa(var,'data.FileMatrix')
 *                 n = var.n;
 *                 m = var.m;
 *                 varargin = {'Dir', fileparts(var.DataDirectory),...
 *                     'BlockSize', var.blocksize}; */
            else
                n = var;
                ip.addRequired('m');
            end

            ip.addParamValue('Dir',KerMor.App.TempDirectory,...
                @(v)isempty(v) || (ischar(v) && exist(v," dir ") == 7));
            ip.addParamValue('BlockSize',KerMor.App.BlockSize,@(v)isposrealscalar(v));
            ip.parse(varargin[:]);
            /*  Latest at here we have an m value */
            if isfield(ip.Results," m ")
                m = ip.Results.m;
            end
            /*  Error checks */
            if ~(round(n) == n && round(m) == m)
                error(" Size arguments n, m must be integer values. ");
            end
            BYTEPERDOUBLE = 8;
            bsize = min(ip.Results.BlockSize*1024^2,n*m*BYTEPERDOUBLE);
            bcols = max(floor(bsize/(BYTEPERDOUBLE*n)),1);
            nb = ceil(m / bcols);

            /*  Determine the DataDirectory (store in each case, even if there is only one block.
             * This is due to the fact that matrix operations can result in a filematrix with
             * more than one block, thus needing to know where the previous matrix should have
             * been located. */
            ddir = ip.Results.Dir;
            if isempty(ddir)
                ddir = KerMor.App.TempDirectory;
            end
            this = this@data.FileData(fullfile(ddir,...
                sprintf(" matrix_%s ",IDGenerator.generateID)));

            this.n= n;
            this.m= m;
            this.blocksize= bsize/1024^2;
            this.bCols= bcols;
            this.nBlocks= nb;
            this.created= false(1,nb);
            hlp = reshape(repmat(1:nb,bcols,1),[],1);
            hlp2 = reshape(repmat(1:bcols,nb,1)^t,[],1);
            this.idx= [hlp(1:m) hlp2(1:m)];

            if nb == 1
                this.loadBlock(1);
            end
            /*  Matrix case: Assign value directly */
            if matrixin
                this.subsasgn(struct(" type ",[" () "]," subs ",[[" : "," : "]]),A);
            end
        }
        function B = spawnWithContent(matrix<double> A) {
            B = data.FileMatrix(A," Dir ",fileparts(this.DataDirectory)," Blocksize ",this.blocksize);
        }
        function [rowmin , rowmax ] = getColBoundingBox() {
            rowmin = Inf*ones(this.n,1);
            rowmax = -rowmin;
            for i=1:this.nBlocks
                A = this.loadBlock(i);
                rowmin = min(rowmin, min(A,[],2));
                rowmax = max(rowmax, max(A,[],2));
            end
        }
        function pos = getBlockPos(nr) {
            pos = (nr-1)*this.bCols + 1 : min(nr*this.bCols,this.m);
        }
        function copy = copyWithNewBlockSize(block_size) {
            copy = data.FileMatrix(this.n,this.m,...
                " Dir ",fileparts(this.DataDirectory)," BlockSize ",block_size);
            for k=1:this.nBlocks
                pos = this.getBlockPos(k);
                copy.subsasgn(struct(" type ",[" () "]," subs ",[[" : ",pos]]),this.loadBlock(k));
            end
        }


        function  relocate(new_root) {

            /*  Extract matrix_xxx folder name */
            [~, mfolder] = fileparts(this.DataDirectory);
            /*  Concatenate with new path and call superclass relocate */
            relocate@data.FileData(this, fullfile(new_root, mfolder));
        }
        function res = transposedTimes(B) {
            if ~isa(B," data.FileMatrix ")
                error(" Operation only applicable to other FileMatrix instances ");
            end
            res = data.FileMatrix(this.m,B.m," Dir ", fileparts(B.DataDirectory),...
                    " BlockSize ", B.blocksize);
            key = struct(" type ",[" () "]," subs ",[]);
            for j = 1:B.nBlocks
                Bb = B.getBlock(j);
                for i=1:this.nBlocks
                    key.subs= [this.getBlockPos(i),B.getBlockPos(j)];
                    res.subsasgn(key,this.getBlock(i)^t*Bb);
                end
            end
        }
        function s = sum(dim) {
            if nargin < 2
                dim = 1;
            end
            if dim == 1
                s = zeros(1,this.m);
                for k=1:this.nBlocks
                    s(this.getBlockPos(k)) = sum(this.loadBlock(k),1);
                end
            elseif dim == 2
                s = zeros(this.n,1);
                for k=1:this.nBlocks
                    s = s + sum(this.loadBlock(k),2);
                end
            else
                error(" Invalid sum dimension: %d ",dim);
            end
        }
        function value = power(expo) {
            if ~isa(this," data.FileMatrix ") && ~isscalar(expo)
                error(" FileMatrix power only defined for scalar values. ");
            end
            /*  Special case for one block */
            if this.nBlocks == 1
                this.ensureCacheLoaded;
                value = this.cachedBlock.^expo;
                return;
            end
            value = data.FileMatrix(this.n,this.m," Dir ",fileparts(this.DataDirectory),...
                " BlockSize ", this.blocksize);
            for k=1:this.nBlocks
                value.saveBlock(k,this.loadBlock(k).^expo);
            end
        }


        function n = numel() {
            n = 1;
        }


        function [n , m ] = size(dim) {
            n = [this.n this.m];
            if nargin == 2
                if dim > 0 && dim < 3
                    n = n(dim);
                else
                    n = 0;
                end
            elseif nargout == 2
                n = this.n;
                m = this.m;
            end
        }
        function varargout = subsref(key) {
            if strcmp(key(1).type," () ")
                if this.nBlocks == 1
                    this.ensureCacheLoaded;
                    [varargout[1:nargout]] = builtin(" subsref ", this.cachedBlock, key);
                    return;
                end
                s = key.subs;
                /*  Default case: row & column adressing */
                if length(s) == 2
                    /*  "Select all" mode */
                    if strcmp(s[2]," : ")
                        s[2] = 1:this.m;
                    end
                    pos = this.idx(s[2],:);
                    blocks = unique(pos(:,1));
                    the_n = this.n;
                    if ~strcmp(s[1]," : ")
                        the_n = length(s[1]);
                    end
                    value = zeros(the_n,size(pos,1));
                    for bidx = 1:length(blocks)
                        b = blocks(bidx);
                        B = this.loadBlock(b);
                        value(:,pos(:,1) == b) = B(s[1],pos(pos(:,1)==b,2));
                    end
                    /*  Linear addressing */
                elseif length(key.subs) == 1
                    error(" Not yet implemented. ");
                else
                    error(" ()-assignment must have one or two subscripts ");
                end
                varargout[1] = value;
            else
                [varargout[1:nargout]] = builtin(" subsref ", this, key);
            end
        }
        function this = subsasgn(key,value) {
            if strcmp(key(1).type," () ")
                if this.nBlocks == 1
                    this.ensureCacheLoaded;
                    this.cachedBlock= builtin(" subsasgn ", this.cachedBlock, key, value);
                    this.cacheDirty= true;
                    return;
                end
                s = key.subs;
                /*  Default case: row & column adressing */
                if length(s) == 2
                    /*  "Select all" mode */
                    if strcmp(s[2]," : ")
                        s[2] = 1:this.m;
                    end
                    pos = this.idx(s[2],:);
                    blocks = unique(pos(:,1));
                    for bidx = 1:length(blocks)
                        b = blocks(bidx);
                        B = this.loadBlock(b);
                        B(s[1],pos(pos(:,1)==b,2)) = value(:,pos(:,1) == b);
                        this.saveBlock(b,B);
                    end
                    /*  Linear addressing */
                elseif length(key.subs) == 1
                    error(" Not yet implemented. ");
                else
                    error(" ()-assignment must have one or two subscripts ");
                end
            else
                builtin(" subsasgn ", this, key, value);
            end
        }
        function diff = minus(B) {
            if all(size(A) == size(B))
                if isa(A," data.FileMatrix ")
                    if isa(B," data.FileMatrix ")
                        error(" Not yet implemented. ");
                        /*                          diff = data.FileMatrix(A.n, A.m, A);
                         *                         key = struct('type',{'()'},'subs',{{':'}});
                         *                         for i=1:A.nBlocks
                         *                             pos = A.getBlockPos(i);
                         *                             key.subs{2} = pos;
                         *                             diff.subsasgn(key,A.loadBlock(i) - B(:,pos));
                         *                         end */
                    else
                        /*  Special case for one block */
                        if A.nBlocks == 1
                            A.ensureCacheLoaded;
                            diff = A.cachedBlock - B;
                            return;
                        end
                        diff = zeros(A.n, A.m);
                        for i=1:A.nBlocks
                            pos = A.getBlockPos(i);
                            diff(:,pos) = A.loadBlock(i) - B(:,pos);
                        end
                    end
                else
                    diff = minus(B, A);
                end
            else
                error(" Invalid matrix dimensions. ");
            end
        }


        function a = mldivide(R) {
            if isa(L," data.FileMatrix ")
                if L.nBlocks == 1
                    L.ensureCacheLoaded;
                    a = L.cachedBlock\R;
                else
                    error(" Not yet implemented ");
                end
            else
                if R.nBlocks == 1
                    R.ensureCacheLoaded;
                    a = L\R.cachedBlock;
                else
                    error(" Not yet implemented ");
                end
            end
        }


        function AB = mtimes(B) {
            if isa(A," data.FileMatrix ")
                /*  FileMatrix * FileMatrix case */
                if isa(B," data.FileMatrix ")
                    if A.nBlocks == 1 && B.nBlocks == 1
                        A.ensureCacheLoaded;
                        B.ensureCacheLoaded;
                        AB = data.FileMatrix(A.cachedBlock * B.cachedBlock);
                        return;
                    else
                        /*  Set blocksize such that the result has the same number of columns per
                         * block as B has */
                        AB = data.FileMatrix(A.n,B.m," Dir ", fileparts(B.DataDirectory),...
                        " BlockSize ", A.blockSizeOf(A.n,B.bCols));
                        key = struct(" type ",[" () "]," subs ",[]);
                        for i=1:B.nBlocks
                            Bb = B.getBlock(i);
                            key.subs= [" : ",B.getBlockPos(i)];
                            for k = 1:A.nBlocks
                                Ab = A.getBlock(k);
                                sum = AB.subsref(key) + Ab*Bb(A.getBlockPos(k),:);
                                AB.subsasgn(key,sum);
                            end
                        end
                    end
                else
                    /*  Special case for one block */
                    if A.nBlocks == 1
                        A.ensureCacheLoaded;
                        AB = A.cachedBlock * B;
                        return;
                    end
                    if isscalar(B)  /*  Matrix*scalar case */

                        AB = mtimes(B,A);
                    else
                        if A.m ~= size(B,1)
                            error(" Matrix dimensions must agree. ");
                        end
                        /*  FileMatrix * vec case */
                        if size(B,2) == 1
                            AB = zeros(A.n,1);
                            for bidx = 1:A.nBlocks
                                /*  Only multiply if nonzero */
                                if A.created(bidx)
                                    ABlock = A.loadBlock(bidx);
                                    AB = AB + ABlock*B(A.getBlockPos(bidx),:);
                                end
                            end
                            /*  FileMatrix * matrix case */
                        else
                            AB = data.FileMatrix(A.n, size(B,2), " Dir ", fileparts(A.DataDirectory),...
                                " BlockSize ", A.blocksize);
                            ABcols = AB.bCols;
                            key = struct(" type ",[" () "]," subs ",[[" : "]]);
                            for i=1:A.nBlocks
                                if A.created(i)
                                    b = A.loadBlock(i);
                                    posm = A.getBlockPos(i);
                                    /*  Split up the additions into the block size of AB */
                                    for j=1:ceil(size(B,2)/ABcols)
                                        key.subs[2] = AB.getBlockPos(j);
                                        /*  This is reasonably fast as the same block in AB is loaded all
                                         * the time. */
                                        AB.subsasgn(key,AB.subsref(key) + b*B(posm,key.subs[2]));
                                    end
                                end
                            end
                        end
                    end
                end
                /*  [vec|mat] * FileMatrix case */
            elseif isa(B," data.FileMatrix ")
                /*  Special case for one block */
                if B.nBlocks == 1
                    B.ensureCacheLoaded;
                    AB = A*B.cachedBlock;
                    return;
                end
                if isscalar(A)  /*  scalar*Matrix */

                    AB = data.FileMatrix(B.n, B.m, " Dir ", fileparts(B.DataDirectory),...
                        " BlockSize ", B.blocksize);
                    key = struct(" type ",[" () "]," subs ",[[" : "]]);
                    for bidx = 1:B.nBlocks
                        /*  Only multiply if nonzero */
                        if B.created(bidx)
                            key.subs[2] = B.getBlockPos(bidx);
                            AB.subsasgn(key,A*B.loadBlock(bidx));
                        end
                    end
                else
                    AB = data.FileMatrix(size(A,1), B.m, " Dir ", fileparts(B.DataDirectory),...
                        " BlockSize ", B.blocksize);
                    key = struct(" type ",[" () "]," subs ",[[" : "]]);
                    /*  If AB has more blocks than B, the resulting matrix is larger than B.
                     * Thus, we need A*b to be as most as big as blocksize, which is why slices of B
                     * in the size of ABs blocks are taken. */
                    if AB.nBlocks > B.nBlocks
                        for i=1:AB.nBlocks
                            key.subs[2] = AB.getBlockPos(i);
                            b = B.subsref(key);
                            AB.subsasgn(key,A*b);
                        end
                    else
                        for i=1:B.nBlocks
                            if B.created(i)
                                key.subs[2] = B.getBlockPos(i);
                                AB.subsasgn(key,A*B.loadBlock(i));
                            end
                        end
                        /*  Return a matrix if the A argument was a transposed vector */
                        if size(A,1) == 1
                            AB = AB.toMemoryMatrix;
                        end
                    end
                end
            end
        }
        function AB = times(B) {
            if isa(A," data.FileMatrix ")
                /*  FileMatrix * FileMatrix case */
                if isa(B," data.FileMatrix ")
                    error(" Not yet implemented. ");
                else
                    /*  Special case for one block */
                    if A.nBlocks == 1
                        A.ensureCacheLoaded;
                        AB = A.cachedBlock .* B;
                        return;
                    elseif isscalar(B)
                        /*  Call matrix multiplication with scalar value */
                        AB = mtimes(B,A);
                    else
                        if  A.n ~= size(B,1) || A.m ~= size(B,2)
                            error(" Matrix dimensions must agree. ");
                        end
                        AB = data.FileMatrix(A.n,A.m," Dir ", fileparts(A.DataDirectory),...
                            " BlockSize ", A.blocksize);

                        for k=1:A.nBlocks
                            AB.saveBlock(k,A.loadBlock(k).*B(:,A.getBlockPos(k)));
                        end
                    end
                end
            else
                /*  Pointwise operation does not depend on which side it is done from. */
                AB = times(B,A);
            end
        }
        function trans = ctranspose() {
            trans = data.FileMatrix(this.m, this.n, " BlockSize ", this.blocksize, ...
                " Dir ", fileparts(this.DataDirectory));
            /*  Special case for only one block */
            if this.nBlocks == 1
                this.ensureCacheLoaded;
                trans = this.cachedBlock^t;
                return;
            end
            /*  Else: multiple blocks */
            if this.InPlaceTranspose
                if this.nBlocks > 1
                    pi = ProcessIndicator(" Creating in-place transposed of %d-block matrix in %s " ,this.nBlocks,...
                        false,this.nBlocks,trans.DataDirectory);
                end
                key = struct(" type ",[" () "]," subs ",[[[]," : "]]);
                for j=1:this.nBlocks
                    B = this.loadBlock(j);
                    key.subs[1] = this.getBlockPos(j);
                    trans.subsasgn(key,B^t);
                    if this.nBlocks > 1
                        pi.step;
                    end
                end
            else
                if this.nBlocks > 1
                    pi = ProcessIndicator(" Creating transposed of %d-block matrix in %s ",2*this.nBlocks,...
                        false,this.nBlocks,trans.DataDirectory);
                end
                /*  Write out blocks in small chunks */
                for j=1:this.nBlocks
                    B = this.loadBlock(j);
                    for k=1:trans.nBlocks
                        pos = trans.getBlockPos(k);
                        chunk = B(pos,:); /* #ok */

                        save(fullfile(trans.DataDirectory,sprintf(" tmp_%d_%d.mat ",j,k)), " chunk ");
                    end
                    if this.nBlocks > 1
                        pi.step;
                    end
                end
                /*  Read in blocks from respective chunks */
                for k=1:trans.nBlocks
                    B = trans.loadBlock(k);
                    for j=1:this.nBlocks
                        file = fullfile(trans.DataDirectory,sprintf(" tmp_%d_%d.mat ",j,k));
                        s = load(file);
                        pos = this.getBlockPos(j);
                        B(pos,:) = s.chunk^t;
                        /*  remove temp file */
                        delete(file);
                    end
                    trans.saveBlock(k,B);
                    if this.nBlocks > 1
                        pi.step;
                    end
                end
            end
            if this.nBlocks > 1
                pi.stop;
            end
        }


        function res = eq(B) {
            res = false;
            if isa(A," data.FileMatrix ") && all(size(A) == size(B))
                if isa(B," data.FileMatrix ")
                    error(" not yet implemented ");
                else
                    for i=1:A.nBlocks
                        pos = A.getBlockPos(i);
                        if A.created(i) || A.cachedNr == i
                            if ~isequal(A.loadBlock(i),B(:,pos)), return; end
                        else
                            if any(any(B(:,pos) ~= 0)), return; end
                        end
                    end
                    res = true;
                end
            elseif isa(B," data.FileMatrix ")
                res = eq(B, A);
            end
        }


        function res = ne(B) {
            res = ~eq(A,B);
        }


        function  horzcat(varargin) {
            error(" not yet implemented ");
        }


        function  vertcat(varargin) {
            error(" not yet implemented ");
        }


        function  delete() {
            if ~isempty(this.created) && ~this.isSaved
                /*  Remove exactly the files for the matrix */
                for k=1:this.nBlocks
                    if this.created(k)
                        f = fullfile(this.DataDirectory, sprintf(" block_%d.mat ",k));
                        if exist(f," file ")
                            delete(f);
                        end
                    end
                end
            end
            /*  Superclass delete removes the folder if empty. */
            delete@data.FileData(this);
        }


        function n = getNumBlocks() {
            n = this.nBlocks;
        }
        function B = getBlock(nr) {
            B = this.loadBlock(nr);
        }
    protected:


        function this = saveobj() {
            this = saveobj@data.FileData(this);
            if this.cacheDirty
                A = this.cachedBlock;/* #ok */

                save([this.DataDirectory filesep sprintf(" block_%d.mat ",this.cachedNr)], " A ");
                this.created(this.cachedNr) = true;
            end
            this.cachedBlock= [];
            this.cachedNr= [];
            this.cacheDirty= false;
        }

    
    private:


        function  ensureCacheLoaded() {
            if this.nBlocks == 1 && isempty(this.cachedBlock)
                this.loadBlock(1);
            end
        }


        function  saveBlock(nr,A) {
            if nr == this.cachedNr
                this.cachedBlock= A;
                this.cacheDirty= true;
            else
                /*  Save actual block that should be saved */
                save([this.DataDirectory filesep sprintf(" block_%d.mat ",nr)], " A ");
            end
            this.updateMinMax(A);
            this.created(nr) = true;
        }
        function A = loadBlock(nr) {
            if nr == this.cachedNr
                A = this.cachedBlock;
            else
                /*  Before loading a new block, save the old one if the cache is dirty! */
                if this.cacheDirty
                    A = this.cachedBlock;
                    save([this.DataDirectory filesep sprintf(" block_%d.mat ",this.cachedNr)], " A ");
                    this.updateMinMax(A);
                    this.created(this.cachedNr) = true;
                    /*  cacheDirty is set "false" at end! */
                end
                if this.created(nr)
                    s = load([this.DataDirectory filesep sprintf(" block_%d.mat ",nr)]);
                    A = s.A;
                else
                    /*  Ensures correct size for block (even if only one with less than bCols
                     * columns) */
                    A = zeros(this.n, length(this.getBlockPos(nr)));
                end
                this.cachedNr= nr;
                this.cachedBlock= A;
                this.cacheDirty= false;
            end
        }


        function  updateMinMax(A) {
            this.MinValue= min(this.MinValue,min(min(A)));
            this.MaxValue= max(this.MaxValue,max(max(A)));
        }

    
    protected: /* ( Static ) */

        static function this = loadobj(this,initfrom) {

            created = false;
            if ~isa(this, " data.FileMatrix ")
                initfrom = this;
                this = data.FileMatrix(initfrom.n,initfrom.m," Dir ",...
                    initfrom.DataDirectory," BlockSize ",initfrom.blocksize);
                created = true;
            end
            if nargin == 2 || created
                this.InPlaceTranspose= initfrom.InPlaceTranspose;
                this.bCols= initfrom.bCols;
                this.nBlocks= initfrom.nBlocks;
                this.n= initfrom.n;
                this.m= initfrom.m;
                this.MinValue= initfrom.MinValue;
                this.MaxValue= initfrom.MaxValue;
                this.idx= initfrom.idx;
                this.created= initfrom.created;
                this.blocksize= initfrom.blocksize;
                if this.nBlocks == 1
                    field = " cachedBlock ";
                    if isfield(initfrom," block1 ")
                        field = " block1 ";
                    end
                    this.cachedBlock= initfrom.(field);
                end
                this = loadobj@data.FileData(this, initfrom);
            else
                this = loadobj@data.FileData(this);
            end
        }
    public: /* ( Static ) */

        static function data.FileMatrixfm = recoverFrom(char directory) {
            if nargin == 0
                directory = KerMor.getDir;
                if ~directory
                    fm = [];
                    return;
                end
            end
            d = dir(directory);
            nf = length(d)-2;
            if nf > 0
                nb = 0;
                for k=1:nf
                    if ~exist(fullfile(directory,sprintf(" block_%d.mat ",k))," file ")
                        break;
                    end
                    nb = nb + 1;
                end
                if nb > 0
                    s = load(fullfile(directory," block_1.mat "));
                    A1 = s.A;
                    [N, nCols] = size(A1);
                    if nb > 1
                        s = load(fullfile(directory,sprintf(" block_%d.mat ",nb)));
                    end
                    M = (nb-1)*nCols + size(s.A,2);
                    fm = data.FileMatrix(N,M," Dir ",fileparts(directory)," BlockSize ",8*nCols*N);
                    fm.updateMinMax(A1); /*  first block */

                    if nb > 1
                        fm.updateMinMax(s.A);
                        for k=2:nb-1
                            s2 = load(fullfile(directory,sprintf(" block_%d.mat ",k)));
                            fm.updateMinMax(s2.A);
                        end
                    end
                    /*  Delete & reassign correct directory */
                    rmdir(fm.DataDirectory);
                    fm.DataDirectory= directory;
                    fm.created(:) = true;
                    fm.isSaved= true;
                else
                    error(" No block_%%.mat files of a FileMatrix found. ");
                end
            else
                error(" No files found in %s ",directory);
            end
        }
        static function bs = blockSizeOf(arg1,integer arg2) {
            if nargin < 2
                n = numel(arg1);
            else
                n = arg1*arg2;
            end
            bs = n*8/1024^2;
        }
        static function res = test_FileMatrix() {
            res = true;
            B = rand(99,100);
            A = data.FileMatrix(99,100," BlockSize ",data.FileMatrix.blockSizeOf(B)/5);
            key = struct(" type ",[" () "]," subs ",[]);
            /*  col-wise setting (fast) */
            for k=1:100
                key.subs= [" : ", k];
                A.subsasgn(key,B(:,k));
            end
            key.subs= [" : ", " : "];
            res = res && all(all(A.subsref(key) == B));
            /*  Row-wise setting */
            for k=1:10
                key.subs= [k," : "];
                A.subsasgn(key,B(k,:));
            end
            key.subs= [" : ", " : "];
            res = res && all(all(A.toMemoryMatrix == B));

            /*  Direct constructor test */
            A2 = data.FileMatrix(B);
            res = res && all(all(A2.subsref(key) == B));

            /*  IsEqual test */
            res = res && A == B;
            B2 = B; B2(1,:) = -B2(1,:);
            res = res && A ~= B2;

            /*  Transpose test */
            A.InPlaceTranspose= false;
            At = A^t;
            res = res && all(all(At.toMemoryMatrix == B^t));
            A.InPlaceTranspose= true;
            At = A^t;
            res = res && all(all(At.toMemoryMatrix == B^t));

            /*  Bounding box test */
            [bm, bM] = Utils.getBoundingBox(B);
            [am, aM] = A.getColBoundingBox;
            res = res && isequal(bm,am) && isequal(bM,aM);
        }


        static function res = test_transposeTimes() {
            [A1,B1] = data.FileMatrix.getTestPair(80,100,5);
            [A2,B2] = data.FileMatrix.getTestPair(80,100,5);
            A3 = A1.transposedTimes(A2);
            res = norm(B1" *B2 - A3.toMemoryMatrix, "fro^t) < sqrt(eps);
        }


        static function res = test_SVD() {
            [A,B] = data.FileMatrix.getTestPair(99,100,5);
            /*  SVD test */
            p = sqrt(eps);
            [u,s,v] = svd(B," econ ");
            [U,S,V] = A.getSVD;
            V = V.toMemoryMatrix;
            res = norm(U*S*V-B," fro ") < p;
            [U5,S5,V5] = A.getSVD(5);
            V5 = V5.toMemoryMatrix;
            res = res && norm(abs(V)-abs(v" ), "fro" ) < p && norm(abs(U)-abs(u), "fro^t) < p &&...
                norm(diag(S)-diag(s)) < p;
            res = res && norm(abs(V5)-abs(v(:,1:5)" ), "fro^t) < p ...
                && norm(abs(U5)-abs(u(:,1:5))," fro ") < p ...
                && norm(diag(S5)-diag(s(1:5,1:5))) < p;

            /*  Exclude test */
            o = general.Orthonormalizer;
            exclu = o.orthonormalize(rand(size(B,1),5));
            [U,~,~] = A.getSVD(10,exclu);
            res = res && norm(exclu^t*U) < 1e-12;
        }


        static function res = test_SaveLoad() {
            [A,B] = data.FileMatrix.getTestPair(100,10000,1);
            /*  Load/save */
            save Atmp A;
            clear A;
            load Atmp;
            res = A == B;
            rmdir(A.DataDirectory," s ");
            delete Atmp.mat;

            /*  Multiblock save/load */
            [A,B] = data.FileMatrix.getTestPair(100,10000,10);

            /*  Load/save */
            save Atmp A;
            clear A;
            load Atmp;
            res = res && A == B;
            rmdir(A.DataDirectory," s ");
            delete Atmp.mat;
        }


        static function res = test_SumPower() {
            [A,B] = data.FileMatrix.getTestPair(100,10000,40);

            /*  Sum test */
            bs = sum(B,2);
            as = sum(A,2);
            res = isequal(sum(B,1),sum(A,1)) && all(abs((as-bs) ./ bs)) < sqrt(eps);

            /*  Power test */
            res = res && A.^2 == B.^2;/* #ok */

        }


        static function res = test_ScalarMult() {
            res = true;
            [A,a] = data.FileMatrix.getTestPair(400,400);
            B = 4*A;
            res = res & isequal(4*a,B);
            C = A*1;
            res = res & isequal(a,C);

            B = 4.*A;
            res = res & isequal(4*a,B);
            C = A.*1;
            res = res & isequal(a,C);

            [A,a] = data.FileMatrix.getTestPair(400,400,4);
            B = 4*A;
            res = res & 4*a == B;
            C = A*1;
            res = res & a == C;

            B = 4.*A;
            res = res & 4*a == B;
            C = A.*1;
            res = res & a == C;
        }
        static function res = test_Times_MTimes() {

            /*  No blocks */
            [A,B] = data.FileMatrix.getTestPair(100,1000);
            AB = A.*B;
            res = isequal(AB,B.*B);
            AB = B.*A;
            res = res & isequal(AB,B.*B);

            /*  With blocks */
            [A,B] = data.FileMatrix.getTestPair(100,1000,4);
            AB = A.*B;
            res = res & AB == B.*B;
            AB = B.*A;
            res = res & AB == B.*B;

            /*  Vector Multiplication test */
            v = rand(size(A,2),1);
            res = res && norm(A*v-B*v) < 1e-10;
            v = rand(size(A,2),100);
            res = res && all(Norm.L2(A*v-B*v) < 1e-10);
            v = rand(1,size(A,1));
            res = res && norm(v*A-v*B) < 1e-10;
            v = rand(100,size(A,1));
            res = res && all(Norm.L2(v*A-v*B) < 1e-10);

            /*  Left / right multiplication test */
            R = rand(1000,200);
            L = rand(2000,50);
            /*  With no blocks */
            [A, Amat] = data.FileMatrix.getTestPair(50,1000);
            LA = L*A;
            AR = A*R;
            res = res && isequal(LA,L*Amat);
            res = res && isequal(AR,Amat*R);
            /*  With blocks */
            [A, Amat] = data.FileMatrix.getTestPair(50,1000,4);
            LA = L*A;
            AR = A*R;
            /*  Left multiplication is exact as whole blocks can be used */
            res = res && LA == L*Amat;
            /*  Right multiplication requires accumulation of values -> rounding errors */
            res = res && norm(AR - Amat*R) < 1e-12;
        }
        static function res = test_2FileMatrix_MTimes() {
            [A,B] = data.FileMatrix.getTestPair(100,1000);
            [C,D] = data.FileMatrix.getTestPair(1000,100);
            AC = A*C;
            res = AC == B*D;
            CA = C*A;
            res = res & CA == D*B;

            /*  With blocks */
            [A,B] = data.FileMatrix.getTestPair(100,1000,4);
            [C,D] = data.FileMatrix.getTestPair(1000,100,5);
            AC = A*C;
            res = res & norm(B*D - AC.toMemoryMatrix," fro ") < 1e-10;
            CA = C*A;
            res = res & norm(D*A - CA.toMemoryMatrix," fro ") < 1e-10;
        }
        static function res = test_SpeedSVDTransp() {
            A = data.FileMatrix.getTestPair(100,1000,5);
            ti = tic;
            [u,s,v] = A.getSVD(5);
            t(1) = toc(ti);

            ti = tic;
            B = A^t;
            t(2) = toc(ti);

            ti = tic;
            [v2,s2,u2] = B.getSVD(5);
            t(3) = toc(ti);
            t(2) = t(2) + t(3);

            fprintf(" Direct SVD time: %g, Transposition time: %g, Transposed SVD time: %g\n ",t);
            res = true;
        }
        static function  test_PartialSVD() {
            d = 100;
            A = data.FileMatrix.getTestPair(d,1000,3);
            B = data.FileMatrix.getTestPair(d,1000,1);
            idx = [1 d-1 floor(rand(1,10)*d)]+1;
            for k = 1:length(idx)
                [u,s] = B.getSVD(100,[],1:idx(k));
                fprintf(" u_nrows=%d, s_nrows=%d\n ",size(u,1),size(s,1));
                [u,s] = A.getSVD(100,[],1:idx(k));
                fprintf(" u_nrows=%d, s_nrows=%d\n ",size(u,1),size(s,1));
            end
        }
    private: /* ( Static ) */

        static function [A , B ] = getTestPair(n,m,nb) {
            if nargin < 3
                nb = 1;
            end
            B = rand(n,m);
            A = data.FileMatrix(n,m," BlockSize ",data.FileMatrix.blockSizeOf(B)/nb);
            /*  Tedious formulation as use of direct operators inside the defining class of
             * overloads will not work.
             * See http://www.mathworks.com/help/matlab/matlab_oop/indexed-reference-and-assignment.html#br09nsm */
            A.subsasgn(struct(" type ",[" () "]," subs ",[[" : ", " : "]]),B);
        }


};
}