Home > utilities > testmextrcolumn.m

# testmextrcolumn

## PURPOSE

TESTMEXTRCOLUMN tests that the behaviour of mextrcolumn() is correct

## SYNOPSIS

function [status] = testmextrcolumn(verbose,tol)

## DESCRIPTION

     TESTMEXTRCOLUMN tests that the behaviour of mextrcolumn() is correct

     RETURN
       status ...  0 failed, 1 all ok
     INPUT
       verbose ... 0 print only what failed, 1 print progress,
                   optional, default 1
       tol ....... tolerance when comparing the results
                   optional, default 1e-8

     See also mextrcolumn, penlabtest, penlabstresstest

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SUBFUNCTIONS

- function [status] = onematrixtest(verbose,tol,testname,A,S1,B,S2,ref)
- function [A,S1,B,Si] = vecmatrixdata(n,nmat,dS1,dSi)
- function [status] = vecmatrixtest(verbose,tol,testname,A,S1,B,Si,ref)

## SOURCE CODE

    0001 function [status] = testmextrcolumn(verbose,tol)
    0002 % TESTMEXTRCOLUMN tests that the behaviour of mextrcolumn() is correct
    0003 %
    0004 % RETURN
    0005 %   status ...  0 failed, 1 all ok
    0006 % INPUT
    0007 %   verbose ... 0 print only what failed, 1 print progress,
    0008 %               optional, default 1
    0009 %   tol ....... tolerance when comparing the results
    0010 %               optional, default 1e-8
    0011 %
    0012 % See also mextrcolumn, penlabtest, penlabstresstest
    0013 %
    0014
    0015 % This file is a part of PENLAB package distributed under GPLv3 license
    0016 % Copyright (c) 2013 by  J. Fiala, M. Kocvara, M. Stingl
    0017 % Last Modified: 5 Dec 2013
    0018
    0019   if (nargin<2)
    0020     tol = 1e-8;
    0021   end
    0022   if (nargin<1)
    0023     verbose = 1;
    0024   end
    0025   status = 1;
    0026
    0027   if (verbose)
    0028     disp(' testing functionality of mextrcolumn()...')
    0029   end
    0030
    0031   % is it possible to call it?
    0032   try
    0033     mextrcolumn(1,sparse(1),1,{sparse(1)});
    0034   catch
    0035     disp(' Error: mextrcolumn - cannot call')
    0036     status = 0;
    0037     return;
    0038   end
    0039
    0040   status = status && onematrixtest(verbose,tol,'n=0, all empty', [],[],[],[],0);
    0041   status = status && onematrixtest(verbose,tol,'n=1, all ones', 1,1,1,1,2);
    0042   status = status && onematrixtest(verbose,tol,'n=1, all nonzeros', 2,3,-4,5,-120*2);
    0043   status = status && onematrixtest(verbose,tol,'n=1, A=0', 0,3,-4,5,0);
    0044   status = status && onematrixtest(verbose,tol,'n=1, S1=0', 2,0,-4,5,0);
    0045   status = status && onematrixtest(verbose,tol,'n=1, B=0', 2,3,0,5,0);
    0046   status = status && onematrixtest(verbose,tol,'n=1, S2=0', 2,3,-4,0,0);
    0047   status = status && onematrixtest(verbose,tol,'n=2, all numbers, A diag', diag([2 3]),[-2.1,3.2;3.2,4.4],[4,6;6,5],[-13,-9;-9,-16]);
    0048   % don't forget that magic() is not symmetric
    0049   status = status && onematrixtest(verbose,tol,'n=2, all numbers, A diag, B=[]', diag([2 3]),[-2.1,3.2;3.2,4.4],magic(2)+magic(2)'+3,sparse(2,2),0);
    0050
    0051   % build special examples with very sparse matrices
    0052   n=5;
    0053   A=diag(1:n);
    0054   B=magic(n)+magic(n)';
    0055   S1=sparse(n,n); S1(2,1)=2;
    0056   S1=S1+S1';
    0057   S2=sparse(n,n); S2(4,4)=9;
    0058   status = status && onematrixtest(verbose,tol,'n=5, very sparse', A,S1,B,S2);
    0059
    0060   n=5;
    0061   A=diag(1:n);
    0062   B=magic(n)+magic(n)';
    0063   S1=sparse(n,n); S1(2,1)=2;
    0064   S1=S1+S1';
    0065   S2=sparse(n,n); S2(2,2)=9;
    0066   status = status && onematrixtest(verbose,tol,'n=5, very sparse2', A,S1,B,S2);
    0067
    0068   n=5;
    0069   A=zeros(n,n);
    0070   A(:)=[1:n^2]-10;
    0071   A=A+A';
    0072   B=magic(n)+magic(n)';
    0073   S1=sparse(n,n); S1(2,1)=2;
    0074   S1=S1+S1';
    0075   S2=sparse(n,n); S2(2,2)=9;
    0076   status = status && onematrixtest(verbose,tol,'n=5, very sparse3', A,S1,B,S2);
    0077
    0078   n=5;
    0079   A=zeros(n,n);
    0080   A(:)=[1:n^2]-10;
    0081   A=A+A';
    0082   B=magic(n)+magic(n)';
    0083   S1=sparse(n,n); S1(2,1)=2;
    0084   S1=S1+S1';
    0085   S2=sparse(diag([5:9]));
    0086   status = status && onematrixtest(verbose,tol,'n=5, sparse, S2 spdiag', A,S1,B,S2);
    0087
    0088   % random of different size
    0089   for n=[5 5 5 10 10 10 100 100 999 1000]
    0090   %for n=[2 5 5]
    0091     A=rand(n,n);
    0092     A=A+A';
    0093     B=rand(n,n);
    0094     B=B+B';
    0095     S1=sprandsym(n,0.2);
    0096     S2=sprandsym(n,0.1);
    0097     testname=sprintf('n=%i all random',n);
    0098     status = status && onematrixtest(verbose,tol,testname, A,S1,B,S2);
    0099   end
    0100
    0101   % test mex on the vector of matrices S2, random data
    0102   %[A,S1,B,Si] = vecmatrixdata(n,nmat,dS1,dSi);
    0103   [A,S1,B,Si] = vecmatrixdata(7,1,0.2,0.9); testname='n=7 * 1 mat';
    0104   status = status && vecmatrixtest(verbose,tol,testname,A,S1,B,Si);
    0105   [A,S1,B,Si] = vecmatrixdata(7,31,0.8,rand(31,1)); testname='n=7 * 31 mat';
    0106   status = status && vecmatrixtest(verbose,tol,testname,A,S1,B,Si);
    0107   [A,S1,B,Si] = vecmatrixdata(8,17,0.5,[17:-1:1]/20); testname='n=8 * 17 mat';
    0108   status = status && vecmatrixtest(verbose,tol,testname,A,S1,B,Si);
    0109   [A,S1,B,Si] = vecmatrixdata(20,10,0.2,ones(10,1)*0.8); testname='n=20 * 10 mat (0.8)';
    0110   status = status && vecmatrixtest(verbose,tol,testname,A,S1,B,Si);
    0111   [A,S1,B,Si] = vecmatrixdata(3,99,0.9,ones(99,1)*0.9); testname='n=3 * 99 mat (0.9)';
    0112   status = status && vecmatrixtest(verbose,tol,testname,A,S1,B,Si);
    0113
    0114 end
    0115
    0116 function [status] = onematrixtest(verbose,tol,testname,A,S1,B,S2,ref)
    0117 % ONEMATRIXTEST - run one test with mextrcolumn(), S2 has just one element
    0118 %   testname - string to print out
    0119 %   A,B,S1,S2 - matrices for the test
    0120 %   ref [optional] - reference result, if not present, computed in Matlab
    0121
    0122   status = 1;
    0123
    0124   if (~issparse(S1))
    0125     S1=sparse(S1);
    0126   end
    0127   if (~issparse(S2))
    0128     S2=sparse(S2);
    0129   end
    0130
    0131   % compute reference result if not provided
    0132   if (nargin==7)
    0133     %ref = trace(A*S1*B*S2);
    0134     ref = trace(A*S1*B*S2) + trace(A*S2*B*S1);
    0135   end
    0136
    0137   try
    0138     result = mextrcolumn(A,S1,B,{S2});
    0139     err = abs(result-ref);
    0140   catch
    0141     err = [];
    0142   end
    0143
    0144   if (isempty(err))
    0145     fprintf(' test %-30s: HARD FAILED\n', testname);
    0146     status = 0;
    0147   elseif (err>tol*abs(ref))
    0148     fprintf(' test %-30s: FAILED  %15.6e %15.6e\n', testname,err,ref);
    0149     status = 0;
    0150   elseif (verbose)
    0151     fprintf(' test %-30s: ok (tol)  %15.6e %15.6e\n', testname,err,ref);
    0152   end
    0153
    0154 end
    0155
    0156 function [A,S1,B,Si] = vecmatrixdata(n,nmat,dS1,dSi)
    0157 % VECMATRIXDATA - generate random dataset for vecmatrixtest
    0158 %   n - dimension of the matrices
    0159 %   nmat - the number of the cell array Si
    0160 %   dS1 - density of S1
    0161 %   dSi - vector of densities (of length nmat) for Si
    0162
    0163   A = rand(n,n);
    0164   A = A+A';
    0165   B = rand(n,n);
    0166   B = B+B';
    0167   S1 = sprandsym(n,dS1);
    0168   Si = cell(nmat,1);
    0169   for i = 1:nmat
    0170     Si{i} = sprandsym(n,dSi(i));
    0171   end
    0172
    0173 end
    0174
    0175 function [status] = vecmatrixtest(verbose,tol,testname,A,S1,B,Si,ref)
    0176 % VECMATRIXTEST - run one test with mextrcolumn() on vector Si
    0177 %   testname - string to print out
    0178 %   A,B,S1,Si{:} - matrices for the test, (Si is a cell array of matrices),
    0179 %      all symmetric, S1, Si{:} must be sparse
    0180 %   ref [optional] - reference result, if not present, computed in Matlab
    0181
    0182   status = 1;
    0183
    0184   % compute reference results if not provided
    0185   if (nargin==7)
    0186     nmat = length(Si);
    0187     ref = zeros(nmat,1);
    0188     for i = 1:nmat
    0189       ref(i) = trace(A*S1*B*Si{i}) + trace(A*Si{i}*B*S1);
    0190     end
    0191   end
    0192   nrm = norm(ref);
    0193
    0194   try
    0195     result = mextrcolumn(A,S1,B,Si);
    0196     err = norm(result-ref);
    0197   catch
    0198     err = [];
    0199   end
    0200
    0201   if (isempty(err))
    0202     fprintf(' test (vec) %-24s: HARD FAILED\n', testname);
    0203     status = 0;
    0204   elseif (err>tol*nrm)
    0205     fprintf(' test (vec) %-24s: FAILED  %15.6e %15.6e\n', testname,err,nrm);
    0206     status = 0;
    0207   elseif (verbose)
    0208     fprintf(' test (vec) %-24s: ok (tol)  %15.6e %15.6e\n', testname,err,nrm);
    0209   end
    0210
    0211 end
    0212