Added the .dia.ref files back in the unit tests.
Michael Baudin [Thu, 7 Aug 2008 12:55:17 +0000 (12:55 +0000)]
36 files changed:
scilab/modules/linear_algebra/tests/unit_tests/balanc.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/balanc.tst
scilab/modules/linear_algebra/tests/unit_tests/bdiag.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/bdiag.tst
scilab/modules/linear_algebra/tests/unit_tests/chol.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/chol.tst
scilab/modules/linear_algebra/tests/unit_tests/companion.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/companion.tst
scilab/modules/linear_algebra/tests/unit_tests/det.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/det.tst
scilab/modules/linear_algebra/tests/unit_tests/gspec.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/gspec.tst
scilab/modules/linear_algebra/tests/unit_tests/hess.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/hess.tst
scilab/modules/linear_algebra/tests/unit_tests/inv.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/inv.tst
scilab/modules/linear_algebra/tests/unit_tests/ldiv.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/ldiv.tst
scilab/modules/linear_algebra/tests/unit_tests/lsq.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/lsq.tst
scilab/modules/linear_algebra/tests/unit_tests/lu.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/lu.tst
scilab/modules/linear_algebra/tests/unit_tests/norm.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/norm.tst
scilab/modules/linear_algebra/tests/unit_tests/qr.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/qr.tst
scilab/modules/linear_algebra/tests/unit_tests/rcond.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/rcond.tst
scilab/modules/linear_algebra/tests/unit_tests/rdiv.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/rdiv.tst
scilab/modules/linear_algebra/tests/unit_tests/schur.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/schur.tst
scilab/modules/linear_algebra/tests/unit_tests/spec.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/spec.tst
scilab/modules/linear_algebra/tests/unit_tests/svd.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/svd.tst

diff --git a/scilab/modules/linear_algebra/tests/unit_tests/balanc.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/balanc.dia.ref
new file mode 100644 (file)
index 0000000..17a8578
--- /dev/null
@@ -0,0 +1,72 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================   balanc     ==============================
+//==========================================================================
+function A=testmat3(a,n)
+//eigen values are given by a dilation of nth roots of 1
+A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+A(1,n)=1/a;
+endfunction
+A=[];[Ab,X]=balanc(A);
+if Ab<>[]|X<>[] then bugmes();quit;end
+// MATRIX
+//Small dimension
+//---------------
+warning('off')
+//Real Case
+A=testmat3(10000,5);
+[Ab,X]=balanc(A);
+if Err(diag(diag(X))-X)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*X)>200*%eps then bugmes();quit;end
+
+//Complex Case
+A=testmat3(10000+0.01*%i,5);
+[Ab,X]=balanc(A);
+if Err(diag(diag(X))-X)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*X)>200*%eps then bugmes();quit;end
+
+//LArge dimension
+A=testmat3(10000,30);
+[Ab,X]=balanc(A);
+if Err(diag(diag(X))-X)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*X)>1000*%eps then bugmes();quit;end
+
+//Complex Case
+A=testmat3(10000+0.01*%i,30);
+[Ab,X]=balanc(A);
+if Err(diag(diag(X))-X)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*X)>1000*%eps then bugmes();quit;end
+
+// PENCILS
+//Small dimension
+//---------------
+//Real Case
+A=testmat3(10000,5);B=testmat3(1000,5);
+[Ab,Bb,X,Y]=balanc(A,B);
+if Err(Bb-inv(X)*B*Y)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*Y)>200*%eps then bugmes();quit;end
+//complex case
+A=testmat3(10000+0.001*%i,5);B=testmat3(1000+100*%i,5);
+[Ab,Bb,X,Y]=balanc(A,B);
+if Err(Bb-inv(X)*B*Y)>200*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*Y)>200*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+//Real Case
+A=testmat3(10000,20);B=testmat3(1000,20);
+[Ab,Bb,X,Y]=balanc(A,B);
+if Err(Bb-inv(X)*B*Y)>1000*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*Y)>1000*%eps then bugmes();quit;end
+//complex case
+A=testmat3(10000+0.001*%i,20);B=testmat3(1000+100*%i,20);
+[Ab,Bb,X,Y]=balanc(A,B);
+if Err(Bb-inv(X)*B*Y)>1000*%eps then bugmes();quit;end
+if Err(Ab-inv(X)*A*Y)>1000*%eps then bugmes();quit;end
+warning('on')
index 21f8965..e49eab8 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/bdiag.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/bdiag.dia.ref
new file mode 100644 (file)
index 0000000..6845db7
--- /dev/null
@@ -0,0 +1,113 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================   bdiag      ==============================
+//==========================================================================
+if bdiag([])<>[] then bugmes();quit;end
+[ab,x]=bdiag([]);
+if ab<>[]|x<>[] then bugmes();quit;end
+[ab,x,bs]=bdiag([]);
+if ab<>[]|x<>[]|bs<>[] then bugmes();quit;end
+if execstr('bdiag([1 2;3 4;5 6])','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+//Real case
+e=1.d-1;
+A=[1 1  2 3 4 5
+   0 1  6 7 8 9
+   0 0  1 e 3 1
+   0 0 -e 1 5 9
+   0 0  0 0 2 e
+   0 0  0 0 0 3];
+X1=[0.5,0.3,0,0.3,0.3,0.2;
+   1,0.6,0.5,0.1,0.7,0.4;
+   0.7,0.1,0.4,0.6,0.1,1;
+   0,0.6,0.2,0.3,0.4,0.5;
+   0.6,0.7,0.5,0.7,0.7,0.5;
+   0.3,0.3,0.4,0.5,0.9,0.6]
+ X1  =
+    0.5    0.3    0.     0.3    0.3    0.2  
+    1.     0.6    0.5    0.1    0.7    0.4  
+    0.7    0.1    0.4    0.6    0.1    1.   
+    0.     0.6    0.2    0.3    0.4    0.5  
+    0.6    0.7    0.5    0.7    0.7    0.5  
+    0.3    0.3    0.4    0.5    0.9    0.6  
+A=inv(X1)*A*X1;
+Ab1=bdiag(A);
+if or(triu(Ab1,-1)<>Ab1) then bugmes();quit;end
+[Ab2,X]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>1d6*%eps then bugmes();quit;end
+[Ab2,X,bs]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>2.d-10 then bugmes();quit;end
+if or(size(bs)<>[3,1]) then bugmes();quit;end
+if sum(bs)<>size(A,1) then bugmes();quit;end
+if or(bs<=0) then bugmes();quit;end
+[Ab2,X,bs]=bdiag(A,1);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>2d-7 then bugmes();quit;end
+if or(size(bs)<>[1,1]) then bugmes();quit;end
+if sum(bs)<>size(A,1) then bugmes();quit;end
+if or(bs<=0) then bugmes();quit;end
+//Complex case
+e=1.d-1;
+A=[1 1  2 3 4 5
+   0 1  6 7 8 9
+   0 0  1 e 3 1
+   0 0 -e 1 5 9
+   0 0  0 0 2 e
+   0 0  0 0 0 3];
+X1=[0.5,0.3,0,0.3,0.3,0.2;
+   1,0.6,0.5,0.1,0.7,0.4;
+   0.7,0.1,0.4,0.6,0.1,1;
+   0,0.6,0.2,0.3,0.4,0.5;
+   0.6,0.7,0.5,0.7,0.7,0.5;
+   0.3,0.3,0.4,0.5,0.9,0.6]+%i*eye(A);
+A=inv(X1)*A*X1;
+Ab1=bdiag(A);
+if or(triu(Ab1)<>Ab1) then bugmes();quit;end
+[Ab2,X]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>1.d-8 then bugmes();quit;end
+[Ab2,X,bs]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>1.d-8 then bugmes();quit;end
+if size(bs,2)<>1 then bugmes();quit;end
+if sum(bs)<>size(A,1) then bugmes();quit;end
+if or(bs<=0) then bugmes();quit;end
+//Large dimension
+//---------------
+//Real case
+A=rand(25,25);
+Ab1=bdiag(A);
+if or(triu(Ab1,-1)<>Ab1) then bugmes();quit;end
+[Ab2,X]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end
+[Ab2,X,bs]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end
+if size(bs,2)<>1 then bugmes();quit;end
+if sum(bs)<>size(A,1) then bugmes();quit;end
+if or(bs<=0) then bugmes();quit;end
+//Complex case
+A=rand(25,25)+%i*rand(25,25);
+Ab1=bdiag(A);
+if or(triu(Ab1)<>Ab1) then bugmes();quit;end
+[Ab2,X]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end
+[Ab2,X,bs]=bdiag(A);
+if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end
+if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end
+if size(bs,2)<>1 then bugmes();quit;end
+if sum(bs)<>size(A,1) then bugmes();quit;end
+if or(bs<=0) then bugmes();quit;end
index b5b78b2..518c235 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/chol.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/chol.dia.ref
new file mode 100644 (file)
index 0000000..b86cf83
--- /dev/null
@@ -0,0 +1,37 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================    chol      ==============================
+//==========================================================================
+//Empty matrix
+if chol([])<>[] then bugmes();quit;end
+if execstr('chol([1 2;3 4])','errcatch')==0 then bugmes();quit;end
+if execstr('chol([1 2;3 4]+%i)','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//REAL
+A=rand(5,5);A=A*A';
+U=chol(A);
+if Err(triu(U)-U)>200*%eps then bugmes();quit;end
+if Err(U'*U-A)>200*%eps then bugmes();quit;end
+//Complex
+A=rand(5,5)+%i*rand(5,5);A=A*A';
+U=chol(A);
+if Err(triu(U)-U)>200*%eps then bugmes();quit;end
+if Err(U'*U-A)>200*%eps then bugmes();quit;end
+//Large dimension
+//REAL
+A=rand(50,50);A=A*A';
+U=chol(A);
+if Err(triu(U)-U)>10*%eps then bugmes();quit;end
+if Err(U'*U-A)>1000*%eps then bugmes();quit;end
+//Complex
+A=rand(5,5)+%i*rand(5,5);A=A*A';
+U=chol(A);
+if Err(triu(U)-U)>10*%eps then bugmes();quit;end
+if Err(U'*U-A)>1000*%eps then bugmes();quit;end
index 5207ee8..2c18be4 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/companion.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/companion.dia.ref
new file mode 100644 (file)
index 0000000..3f57ad5
--- /dev/null
@@ -0,0 +1,55 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// companion.tst --
+//   Test "companion" for real and complex polynomials
+//   with linear, quadratic and cubic polynomials.
+//
+// Linear real polynomial
+p=1+2*%s;
+computed=companion(p);
+expected=[-1/2];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Quadratic real polynomial
+p=1+2*%s+3*%s^2;
+computed=companion(p);
+expected=[-2/3 , -1/3;1 , 0];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Cubic real polynomial
+p=1+2*%s+3*%s^2+4*%s^3;
+computed=companion(p);
+expected=[-3/4 , -2/4 , -1/4; 1 , 0 , 0 ; 0 , 1, 0];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Linear complex polynomial
+p=1+%i+2*%s;
+computed=companion(p);
+expected=[-(1+%i)/2];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Quadratic complex polynomial
+p=1+%i+2*%s+3*%s^2;
+computed=companion(p);
+expected=[-2/3 , -(1+%i)/3;1 , 0];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Cubic complex polynomial
+p=1+%i+2*%s+3*%s^2+4*%s^3;
+computed=companion(p);
+expected=[-3/4 , -2/4 , -(1+%i)/4; 1 , 0 , 0 ; 0 , 1, 0];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Vector of linear polynomials
+p1=1+2*%s;
+p2=1+%i+2*%s;
+vector = [p1 p2];
+computed=companion(vector);
+expected=[-1/2 0;0 -(1+%i)/2];
+if abs(expected-computed)>%eps then bugmes();quit;end
+// Vector of quadratic/cubic real/complex polynomials
+p1=1+2*%s+3*%s^2;
+p2=1+%i+2*%s+3*%s^2+4*%s^3;
+vector = [p1 p2];
+computed=companion(vector);
+expected=[-2/3 -1/3 0 0 0;1 0 0 0 0; 0 0 -3/4 -2/4 -(1+%i)/4;0 0 1 0 0;0 0 0 1 0];
+if abs(expected-computed)>%eps then bugmes();quit;end
index 1d7d836..68d58cc 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 //
 // companion.tst --
 //   Test "companion" for real and complex polynomials
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/det.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/det.dia.ref
new file mode 100644 (file)
index 0000000..c98bd71
--- /dev/null
@@ -0,0 +1,36 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================   det        ==============================
+//==========================================================================
+if execstr('det([1 2;3 4;5 6])','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//Real
+A=[1 1; 1 2];
+if Err(det(A)-1)>10*%eps then bugmes();quit;end
+[e,m]=det(A);
+if e<>0 |Err(m-1)>10*%eps then bugmes();quit;end
+//Complex
+A=A+%i;
+if Err(det(A)-1-%i)>10*%eps then bugmes();quit;end
+[e,m]=det(A);
+if e<>0 |Err(m-1-%i)>10*%eps then bugmes();quit;end
+//Large dimension
+//Real
+v=rand(1,21);
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+if Err(det(A)-prod(v))>400000*%eps then bugmes();quit;end
+[e,m]=det(A);
+if Err(m*(10^e)-prod(v))>400000*%eps then bugmes();quit;end
+//Complex
+v=(v+rand(v)*%i)/2;
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+if Err(det(A)-prod(v))>10000*%eps then bugmes();quit;end
+[e,m]=det(A);
+if Err(m*(10^e)-prod(v))>10000*%eps then bugmes();quit;end
index 0e90b6a..c8f9c39 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/gspec.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/gspec.dia.ref
new file mode 100644 (file)
index 0000000..1e774d4
--- /dev/null
@@ -0,0 +1,83 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+rand('normal')
+//define tools
+function A=testmat1(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+       A(1,n)=1/a;A(n,1)=a
+endfunction
+function r=Err(x)
+       r=norm(x,1)
+endfunction
+//==========================================================================
+//==============================    gspec     ==============================
+//==========================================================================
+//Empty matrix
+S=spec([],[]);
+if S<>[] then bugmes();quit;end
+[al,be]=spec([],[]);
+if al<>[]|be<>[] then bugmes();quit;end
+[al,be,Z]=spec([],[]);
+if al<>[]|be<>[]|Z<>[] then bugmes();quit;end
+[al,be,Z,Q]=spec([],[]);
+if al<>[]|be<>[]|Z<>[]|Q<>[] then bugmes();quit;end
+//Matrix with Inf or Nan (test de la detection d'erreur
+if execstr('spec([%inf 1;2 3],[1 2;3 4])','errcatch')==0 then bugmes();quit;end
+if execstr('spec([1 2;3 4],[1 %nan;2 3])','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+//Real Case
+A=testmat1(3,5);E=testmat1(-2,5);
+S=spec(A,E);
+[Sa,Se]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+[Sa,Se,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>200*%eps then bugmes();quit;end
+[Sa,Se,Q,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>200*%eps then bugmes();quit;end
+if Err(Q'*A-diag(S)*Q'*E)>200*%eps then bugmes();quit;end
+//Complex Case
+A=testmat1(3-%i,5);E=testmat1(-2+0.1*%i,5);
+S=spec(A,E);
+[Sa,Se]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+[Sa,Se,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>200*%eps then bugmes();quit;end
+[Sa,Se,Q,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>200*%eps then bugmes();quit;end
+if Err(Q'*A-diag(S)*Q'*E)>200*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+//Real Case
+A=testmat1(3,30);E=testmat1(-2,30);
+S=spec(A,E);
+[Sa,Se]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+[Sa,Se,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>1000*%eps then bugmes();quit;end
+[Sa,Se,Q,Z]=spec(A,E);
+if Err(S-Sa./Se)>10*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>1000*%eps then bugmes();quit;end
+if Err(Q'*A-diag(S)*Q'*E)>1000*%eps then bugmes();quit;end
+//Complex Case
+A=testmat1(3-%i,30);E=testmat1(-2+0.1*%i,30);
+S=spec(A,E);
+[Sa,Se]=spec(A,E);
+if Err(S-Sa./Se)>20*%eps then bugmes();quit;end
+[Sa,Se,Z]=spec(A,E);
+if Err(S-Sa./Se)>20*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>1000*%eps then bugmes();quit;end
+[Sa,Se,Q,Z]=spec(A,E);
+if Err(S-Sa./Se)>20*%eps then bugmes();quit;end
+if Err(A*Z-E*Z*diag(S))>1000*%eps then bugmes();quit;end
+if Err(Q'*A-diag(S)*Q'*E)>1000*%eps then bugmes();quit;end
index 8fe4602..dcf6a75 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 rand('normal')
 
 //define tools
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/hess.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/hess.dia.ref
new file mode 100644 (file)
index 0000000..ceff6d2
--- /dev/null
@@ -0,0 +1,48 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================     hess     ==============================
+//==========================================================================
+//Empty matrix
+if hess([])<>[] then bugmes();quit;end
+[U,H]=hess([]);
+if U<>[]|H<>[] then bugmes();quit;end
+if execstr('hess(rand(2,5))','errcatch')==0 then bugmes();quit;end
+if execstr('[U,H]=hess(rand(2,5))','errcatch')==0 then bugmes();quit;end
+if execstr('hess(rand(2,5)+%i)','errcatch')==0 then bugmes();quit;end
+if execstr('[U,H]=hess(rand(2,5)+%i)','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//Real case
+A=rand(5,5);
+H=hess(A);
+[U,H1]=hess(A);
+if Err(H-H1)>200*%eps then bugmes();quit;end
+if Err(U'*U-eye()) >200*%eps then bugmes();quit;end
+if Err(U'*A*U-H1)  >200*%eps then bugmes();quit;end
+//complex case
+A=rand(5,5)+%i*rand(5,5);
+H=hess(A);
+[U,H1]=hess(A);
+if Err(H-H1)>200*%eps then bugmes();quit;end
+if Err(U'*U-eye()) >200*%eps then bugmes();quit;end
+if Err(U'*A*U-H1)  >200*%eps then bugmes();quit;end
+//Large dimension
+A=rand(20,20);
+H=hess(A);
+[U,H1]=hess(A);
+if Err(H-H1)>200*%eps then bugmes();quit;end
+if Err(U'*U-eye()) >1000*%eps then bugmes();quit;end
+if Err(U'*A*U-H1)  >1000*%eps then bugmes();quit;end
+//complex case
+A=rand(20,20)+%i*rand(20,20);
+H=hess(A);
+[U,H1]=hess(A);
+if Err(H-H1)>1000*%eps then bugmes();quit;end
+if Err(U'*U-eye()) >1000*%eps then bugmes();quit;end
+if Err(U'*A*U-H1)  >1000*%eps then bugmes();quit;end
index f95b346..8fce2bc 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/inv.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/inv.dia.ref
new file mode 100644 (file)
index 0000000..aaad5e3
--- /dev/null
@@ -0,0 +1,56 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//define tools
+function A=testmat1(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+       A(1,n)=1/a;A(n,1)=a
+endfunction
+function r=Err(x)
+       r=norm(x,1)
+endfunction
+rand('normal')
+//==========================================================================
+//==============================     inv      ==============================
+//==========================================================================
+//Empty matrix
+A=[];
+if inv(A)<>[] then bugmes();quit;end
+//Singular matrix
+if execstr('inv([0 0;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('inv([0 0;%i 3])','errcatch')==0 then bugmes();quit;end
+//Rectangular matrix
+if execstr('inv(rand(2,3))','errcatch')==0 then bugmes();quit;end
+if execstr('inv(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+//Unsymetric
+A=testmat1(3,5);Ac=testmat1(3+%i,5);
+//Real Case
+if Err(A*inv(A)-eye(A)) >200*%eps then bugmes();quit;end
+//Complex Case
+if Err(Ac*inv(Ac)-eye(A)) >200*%eps then bugmes();quit;end
+//Symetric
+A=A*A';Ac=Ac*Ac';
+//Real Case
+if Err(A*inv(A)-eye(A)) >1000*%eps then bugmes();quit;end
+//Complex Case
+if Err(Ac*inv(Ac)-eye(A)) >1000*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+//Unsymetric
+A=testmat1(3,50);Ac=testmat1(3+%i,50);
+//Real Case
+if Err(A*inv(A)-eye(A)) >1000*%eps then bugmes();quit;end
+//Complex Case
+if Err(Ac*inv(Ac)-eye(A)) >2000*%eps then bugmes();quit;end
+//Symetric
+A=A*A';Ac=Ac*Ac';
+//Real Case
+if Err(A*inv(A)-eye(A)) >1.d-10 then bugmes();quit;end
+//Complex Case
+if Err(Ac*inv(Ac)-eye(A)) >4.d-10 then bugmes();quit;end
index 0c78d9f..df88e0f 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 //define tools
 function A=testmat1(a,n)
        //eigen values are given by a dilation of nth roots of 1
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/ldiv.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/ldiv.dia.ref
new file mode 100644 (file)
index 0000000..8b036e0
--- /dev/null
@@ -0,0 +1,121 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//define tools
+function A=testmat1(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+       A(1,n)=1/a;A(n,1)=a
+endfunction
+function r=Err(x)
+       r=norm(x,1)
+endfunction
+rand('normal')
+//==========================================================================
+//==============================    \         ==============================
+//==========================================================================
+function x=ldiv(A,B),x=A\B,endfunction
+Warning : redefining function: ldiv                    . Use funcprot(0) to avoid this message
+
+//scalar division
+//-----Square matrix-----
+n=5;  A=rand(n,n);b=rand(2,n+1);
+if ldiv(A,[]) <>[] then bugmes();quit;end
+if execstr('ldiv(A,B)','errcatch')==0 then bugmes();quit;end
+//Small dimensions real
+n=5;
+b=rand(n,2);A=rand(n,n);
+x=ldiv(A,b);
+if Err(A*x-b)>200*%eps then bugmes();quit;end
+//Small dimensions complex
+b=rand(n,2)+%i;A=rand(n,n);
+x=ldiv(A,b);
+if Err(A*x-b)>500*%eps then bugmes();quit;end
+b=rand(n,2);A=rand(n,n)+%i;
+x=ldiv(A,b);
+if Err(A*x-b)>200*%eps then bugmes();quit;end
+b=rand(n,2)+%i;A=rand(n,n)+%i;
+x=ldiv(A,b);
+if Err(A*x-b)>200*%eps then bugmes();quit;end
+//Large dimensions real
+n=50;
+b=rand(n,2);A=rand(n,n);
+x=ldiv(A,b);
+if Err(A*x-b)>10^6*%eps then bugmes();quit;end
+//Small dimensions complex
+b=rand(n,2)+%i;A=rand(n,n);
+x=ldiv(A,b);
+if Err(A*x-b)>50000*%eps then bugmes();quit;end
+b=rand(n,2);A=rand(n,n)+%i;
+x=ldiv(A,b);
+if Err(A*x-b)>50000*%eps then bugmes();quit;end
+b=rand(n,2)+%i;A=rand(n,n)+%i;
+x=ldiv(A,b);
+if Err(A*x-b)>50000*%eps then bugmes();quit;end
+//-----Rectangular matrix-----
+n=5;m=3; A=rand(m,n);b=rand(n+1,2);
+if ldiv(A,[]) <>[] then bugmes();quit;end
+if execstr('ldiv(A,b)','errcatch')==0 then bugmes();quit;end
+//Small dimensions real
+n=5;m=3;
+b=rand(m,2);A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>200*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(m,2);A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+//Small dimensions complex
+n=5;m=3;
+b=rand(m,2)+%i;A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+n=5;m=3;
+b=rand(m,2);A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+b=rand(m,2)+%i;A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(m,2)+%i;A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(m,2);A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(m,2)+%i;A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>1000*%eps then bugmes();quit;end
+//LArge dimension real
+n=40;m=20;
+b=rand(m,2);A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2);A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+//Large dimensions complex
+b=rand(m,2)+%i;A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2);A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2)+%i;A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2)+%i;A=rand(m,n);
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2);A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
+b=rand(m,2)+%i;A=rand(m,n)+%i;
+x=ldiv(A,b);
+if Err(A'*A*x-A'*b)>10000*%eps then bugmes();quit;end
index ab2052d..d05658e 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 //define tools
 function A=testmat1(a,n)
        //eigen values are given by a dilation of nth roots of 1
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/lsq.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/lsq.dia.ref
new file mode 100644 (file)
index 0000000..5fefb8a
--- /dev/null
@@ -0,0 +1,96 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================     lsq      ==============================
+//==========================================================================
+//Empty matrix
+if lsq([],[])<>[] then bugmes();quit;end
+if execstr('lsq([],1)','errcatch')==0 then bugmes();quit;end
+if execstr('lsq(1,[])','errcatch')==0 then bugmes();quit;end
+if execstr('lsq(rand(3,2),rand(2,1))','errcatch')==0 then bugmes();quit;end
+//Small dimensions
+//Real full rank fat
+A=rand(3,5);b=rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+//Real rank deficient fat
+A=[1 2 3;1 2 3];b=[4;5];
+X=lsq(A,b);
+if Err(A'*A*X-A'*b)> 200*%eps then bugmes();quit;end
+//Real  tall
+A=[1 2;4 2;0 1];b=[1;1;1];
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:2,:))*b1(1:2);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+//Complex full rank fat
+A=rand(3,5)+%i*rand(3,5);b=rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+A=rand(3,5);b=rand(3,2)+%i*rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+A=rand(3,5)+%i*rand(3,5);b=rand(3,2)+%i*rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+//Complex  rank deficient fat
+A=[1 2 3;1 2 3]+%i;b=[4;5];
+X=lsq(A,b);
+A=[1 2 3;1 2 3];b=[4;5]+%i;
+X=lsq(A,b);
+if Err(A'*A*X-A'*b)>200*%eps then bugmes();quit;end
+if Err(A'*A*X-A'*b)>200*%eps then bugmes();quit;end
+A=[1 2 3;1 2 3]+%i;b=[4;5]+%i;
+X=lsq(A,b);
+if Err(A'*A*X-A'*b)>1000*%eps then bugmes();quit;end
+//Complex  full rank tall
+A=[1 2;4 2;0 1]+%i;b=[1;1;1];
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:2,:))*b1(1:2);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+A=[1 2;4 2;0 1];b=[1;1;1]+%i;
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:2,:))*b1(1:2);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+A=[1 2;4 2;0 1]+%i;b=[1;1;1]+%i;
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:2,:))*b1(1:2);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+//LArge dimension
+//Real full rank fat
+A=rand(3,50);b=rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+//Real full rank tall
+A=rand(50,3);b=rand(50,2);
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:3,:))*b1(1:3,:);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+//Complex full rank fat
+A=rand(3,50)+%i;b=rand(3,2);
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+A=rand(3,50);b=rand(3,2)+%i;
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+A=rand(3,50);b=rand(3,2)+%i;
+X=lsq(A,b);
+if Err(A*X-b)>200*%eps then bugmes();quit;end
+//Complex full rank tall
+A=rand(50,3)+%i;b=rand(50,2);
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:3,:))*b1(1:3,:);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+A=rand(50,3);b=rand(50,2)+%i;
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:3,:))*b1(1:3,:);
+if Err(X-X1)>200*%eps then bugmes();quit;end
+A=rand(50,3)+%i;b=rand(50,2)+%i;
+X=lsq(A,b);
+[Q,R]=qr(A);b1=Q'*b;X1=inv(R(1:3,:))*b1(1:3,:);
+if Err(X-X1)>200*%eps then bugmes();quit;end
index 462341f..8123a61 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/lu.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/lu.dia.ref
new file mode 100644 (file)
index 0000000..23a6142
--- /dev/null
@@ -0,0 +1,99 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================      lu      ==============================
+//==========================================================================
+//Empty matrix
+A=[];
+[L,U]=lu(A);
+if L<>[]|U<>[] then bugmes();quit;end
+[L,U,E]=lu(A);
+if L<>[]|U<>[]|E<>[] then bugmes();quit;end
+//Non full rank
+A=rand(5,2);A=A*A';;Ac=rand(5,2)+%i*rand(5,2);Ac=Ac*Ac';
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >200*%eps then bugmes();quit;end
+//Small dimension
+//---------------
+//Square
+A=rand(5,5);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >200*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >200*%eps then bugmes();quit;end
+//Fat
+A=rand(3,5);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >200*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >200*%eps then bugmes();quit;end
+//Tall
+A=rand(5,3);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >200*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >200*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >200*%eps then bugmes();quit;end
+//large dimension
+//---------------
+//Square
+A=rand(50,50);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >1000*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >1000*%eps then bugmes();quit;end
+//Fat
+A=rand(30,50);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >1000*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >1000*%eps then bugmes();quit;end
+//Tall
+A=rand(50,30);Ac=A+%i*rand(A);
+//Real case
+[L,U]=lu(A);
+if Err(L*U-A) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(A);
+if Err(L*U-E*A) >1000*%eps then bugmes();quit;end
+//Complex case
+[L,U]=lu(Ac);
+if Err(L*U-Ac) >1000*%eps then bugmes();quit;end
+[L,U,E]=lu(Ac);
+if Err(L*U-E*Ac) >1000*%eps then bugmes();quit;end
index c8e9c99..2bed8cd 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/norm.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/norm.dia.ref
new file mode 100644 (file)
index 0000000..e0af04a
--- /dev/null
@@ -0,0 +1,61 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+eps=100*%eps;
+// inf
+if norm([1,2,3,-1,-2,-3],0)<>%inf then bugmes();quit;end
+if ~isnan(norm([1,2,3,-1,-2,-3],%nan)) then bugmes();quit;end
+if norm([])<>0 then bugmes();quit;end
+// vector
+x=[1,2,3,-4];
+if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
+if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
+if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
+p=0.5;
+if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
+p=2.5;
+if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
+if abs(norm(x,'inf') -maxi(abs(x))) > eps then bugmes();quit;end
+if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
+if abs(norm(x,'fro') -norm(x,2)) > eps then bugmes();quit;end
+// complex
+x=x+%i*x;
+if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
+if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
+if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
+p=0.5;
+// 100*%eps is needed for linux
+if abs(norm(x,p) - maxi(abs(x))*sum((abs(x)/maxi(abs(x)))^p)^(1/p))> 100*%eps then bugmes();quit;end
+p=2.5;
+if abs(norm(x,p) - maxi(abs(x))*sum((abs(x)/maxi(abs(x)))^p)^(1/p))> 100*%eps then bugmes();quit;end
+if abs(norm(x,'inf') -maxi(abs(x)))> eps then bugmes();quit;end
+if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
+if abs(norm(x,'fro') -norm(x,2))> eps then bugmes();quit;end
+// scalar
+x=[1.23];
+if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
+if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
+if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
+p=0.5;
+if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
+p=2.5;
+if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
+if abs(norm(x,'inf') -maxi(abs(x))) > eps then bugmes();quit;end
+if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
+if abs(norm(x,'fro') -norm(x,2)) > eps then bugmes();quit;end
+// Matrices
+a=rand(10,10,'g');
+if abs(norm(a,1) - maxi(sum(abs(a),'r'))) > eps then bugmes();quit;end
+if abs(norm(a,'inf') - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
+if abs(norm(a,%inf) - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
+if abs(norm(a,2) - maxi(svd(a))) > eps then bugmes();quit;end
+if abs(norm(a,'fro') - norm(matrix(a,1,size(a,'*')),2)) > eps then bugmes();quit;end
+a=a+%i*a;
+if abs(norm(a,1) - maxi(sum(abs(a),'r'))) > eps then bugmes();quit;end
+if abs(norm(a,'inf') - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
+if abs(norm(a,%inf) - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
+if abs(norm(a,2) - maxi(svd(a))) > eps then bugmes();quit;end
+if abs(norm(a,'fro') - norm(matrix(a,1,size(a,'*')),2)) > eps then bugmes();quit;end
index 86a36e3..1d90934 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 eps=100*%eps;
 // inf 
 if norm([1,2,3,-1,-2,-3],0)<>%inf then pause,end 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/qr.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/qr.dia.ref
new file mode 100644 (file)
index 0000000..925bf47
--- /dev/null
@@ -0,0 +1,143 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================      qr      ==============================
+//==========================================================================
+//Empty matrix
+e=[];
+if qr(e)<>[] then bugmes();quit;end
+if qr(e,"e")<>[] then bugmes();quit;end
+[Q,R]=qr(e);
+if Q<>[]|R<>[] then bugmes();quit;end
+[Q,R]=qr(e,"e");
+if Q<>[]|R<>[] then bugmes();quit;end
+[Q,R,x]=qr(e);
+if Q<>[]|R<>[]|x<>[] then bugmes();quit;end
+[Q,R,x]=qr(e,"e");
+if Q<>[]|R<>[]|x<>[] then bugmes();quit;end
+//Small dimension
+//---------------
+A=rand(3,2);Ac=A+rand(A)*%i;
+//Real Case
+Q=qr(A);
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+Q=qr(A,"e");
+if Err(Q'*Q-eye())> 200*%eps then bugmes();quit;end
+[Q,R]=qr(A);
+if Err(Q*R-A)> 200*%eps then bugmes();quit;end
+[Q,R]=qr(A,"e");
+if Err(Q*R-A)> 200*%eps then bugmes();quit;end
+if Err(Q'*Q-eye())> 200*%eps then bugmes();quit;end
+Q=qr(A');
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+Q=qr(A',"e");
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+[Q,R]=qr(A');
+if Err(Q*R-A')> 200*%eps then bugmes();quit;end
+[Q,R]=qr(A',"e");
+if Err(Q*R-A')> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(A);
+if Err(Q*R*x'-A)> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(A,"e");
+if Err(Q*R*x'-A)> 200*%eps then bugmes();quit;end
+//Complex case
+Q=qr(Ac);
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+Q=qr(Ac,"e");
+if Err(Q'*Q-eye())> 200*%eps then bugmes();quit;end
+[Q,R]=qr(Ac);
+if Err(Q*R-Ac)> 200*%eps then bugmes();quit;end
+[Q,R]=qr(Ac,"e");
+if Err(Q*R-Ac)> 200*%eps then bugmes();quit;end
+if Err(Q'*Q-eye())> 200*%eps then bugmes();quit;end
+Q=qr(Ac');
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+Q=qr(Ac',"e");
+if Err(Q*Q'-eye())> 200*%eps then bugmes();quit;end
+[Q,R]=qr(Ac');
+if Err(Q*R-Ac')> 200*%eps then bugmes();quit;end
+[Q,R]=qr(Ac',"e");
+if Err(Q*R-Ac')> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac);
+if Err(Q*R-Ac*x)> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac,"e");
+if Err(Q*R-Ac*x)> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac');
+if Err(Q*R-Ac'*x)> 200*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac',"e");
+if Err(Q*R-Ac'*x)> 200*%eps then bugmes();quit;end
+//Rank detection (obsolete)
+[Q,R,rk,x]=qr(A);
+if Err(Q*R*x'-A)> 200*%eps | rk<>2  then bugmes();quit;end
+[Q,R,rk,x]=qr(A,1.d-8);
+if Err(Q*R*x'-A)> 200*%eps | rk<>2  then bugmes();quit;end
+[Q,R,rk,x]=qr(Ac);
+if Err(Q*R*x'-Ac)> 200*%eps | rk<>2  then bugmes();quit;end
+[Q,R,rk,x]=qr(Ac,1.d-8);
+if Err(Q*R*x'-Ac)> 200*%eps | rk<>2  then bugmes();quit;end
+//Large dimension
+//---------------
+A=rand(150,60);Ac=A+rand(A)*%i;
+//Real Case
+Q=qr(A);
+if Err(Q*Q'-eye())> 1000*%eps then bugmes();quit;end
+Q=qr(A,"e");
+if Err(Q'*Q-eye())> 1000*%eps then bugmes();quit;end
+[Q,R]=qr(A);
+if Err(Q*R-A)> 1000*%eps then bugmes();quit;end
+[Q,R]=qr(A,"e");
+if Err(Q*R-A)> 1000*%eps then bugmes();quit;end
+if Err(Q'*Q-eye())> 1000*%eps then bugmes();quit;end
+Q=qr(A');
+if Err(Q*Q'-eye())> 1000*%eps then bugmes();quit;end
+Q=qr(A',"e");
+if Err(Q*Q'-eye())> 1000*%eps then bugmes();quit;end
+[Q,R]=qr(A');
+if Err(Q*R-A')> 1000*%eps then bugmes();quit;end
+[Q,R]=qr(A',"e");
+if Err(Q*R-A')> 1000*%eps then bugmes();quit;end
+[Q,R,x]=qr(A);
+if Err(Q*R*x'-A)> 1000*%eps then bugmes();quit;end
+[Q,R,x]=qr(A,"e");
+if Err(Q*R*x'-A)> 1000*%eps then bugmes();quit;end
+//Complex case
+Q=qr(Ac);
+if Err(Q*Q'-eye())> 2000*%eps then bugmes();quit;end
+Q=qr(Ac,"e");
+if Err(Q'*Q-eye())> 2000*%eps then bugmes();quit;end
+[Q,R]=qr(Ac);
+if Err(Q*R-Ac)> 2000*%eps then bugmes();quit;end
+[Q,R]=qr(Ac,"e");
+if Err(Q*R-Ac)> 2000*%eps then bugmes();quit;end
+if Err(Q'*Q-eye())> 2000*%eps then bugmes();quit;end
+Q=qr(Ac');
+if Err(Q*Q'-eye())> 2000*%eps then bugmes();quit;end
+Q=qr(Ac',"e");
+if Err(Q*Q'-eye())> 2000*%eps then bugmes();quit;end
+[Q,R]=qr(Ac');
+if Err(Q*R-Ac')> 2000*%eps then bugmes();quit;end
+[Q,R]=qr(Ac',"e");
+if Err(Q*R-Ac')> 2000*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac);
+if Err(Q*R-Ac*x)> 2000*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac,"e");
+if Err(Q*R-Ac*x)> 2000*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac');
+if Err(Q*R-Ac'*x)> 2000*%eps then bugmes();quit;end
+[Q,R,x]=qr(Ac',"e");
+if Err(Q*R-Ac'*x)> 2000*%eps then bugmes();quit;end
+//Rank detection (obsolete)
+[Q,R,rk,x]=qr(A);
+if Err(Q*R*x'-A)> 2000*%eps | rk<>60  then bugmes();quit;end
+[Q,R,rk,x]=qr(A,1.d-8);
+if Err(Q*R*x'-A)> 2000*%eps | rk<>60  then bugmes();quit;end
+[Q,R,rk,x]=qr(Ac);
+if Err(Q*R*x'-Ac)> 2000*%eps | rk<>60  then bugmes();quit;end
+[Q,R,rk,x]=qr(Ac,1.d-8);
+if Err(Q*R*x'-Ac)> 2000*%eps | rk<>60  then bugmes();quit;end
index 76abbd3..8bdf1c8 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 //==========================================================================
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/rcond.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/rcond.dia.ref
new file mode 100644 (file)
index 0000000..28de8c3
--- /dev/null
@@ -0,0 +1,29 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================    rcond     ==============================
+//==========================================================================
+//Empty matrix
+A=[];
+if rcond(A)<>[] then bugmes();quit;end
+//Rectangular matrix
+if execstr('rcond(rand(2,3))','errcatch')==0 then bugmes();quit;end
+if execstr('rcond(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+//Real Case
+if Err(rcond(eye(5,5))-1)>10*%eps then bugmes();quit;end
+//Complex Case
+if  Err(rcond(eye(5,5)*(1+%i))-1)>10*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+//Real Case
+if Err(rcond(eye(50,50))-1)>10*%eps then bugmes();quit;end
+//Complex Case
+if  Err(rcond(eye(50,50)*(1+%i))-1)>10*%eps then bugmes();quit;end
index c4ea070..ca2b610 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/rdiv.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/rdiv.dia.ref
new file mode 100644 (file)
index 0000000..ee39fc4
--- /dev/null
@@ -0,0 +1,118 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================    /         ==============================
+//==========================================================================
+//function x=rdiv(A,B),x=slash(A,B),endfunction
+function x=rdiv(A,B),x=A/B,endfunction
+//scalar division
+//-----Square matrix-----
+n=5;  A=rand(n,n);b=rand(2,n+1);
+if rdiv([],A) <>[] then bugmes();quit;end
+if execstr('rdiv(b,A)','errcatch')==0 then bugmes();quit;end
+//Small dimensions real
+n=5;
+b=rand(2,n);A=rand(n,n);
+x=rdiv(b,A);
+if Err(x*A-b)>200*%eps then bugmes();quit;end
+//Small dimensions complex
+b=rand(2,n)+%i;A=rand(n,n);
+x=rdiv(b,A);
+if Err(x*A-b)>500*%eps then bugmes();quit;end
+b=rand(2,n);A=rand(n,n)+%i;
+x=rdiv(b,A);
+if Err(x*A-b)>500*%eps then bugmes();quit;end
+b=rand(2,n)+%i;A=rand(n,n)+%i;
+x=rdiv(b,A);
+if Err(x*A-b)>500*%eps then bugmes();quit;end
+//Large dimensions real
+n=50;
+b=rand(2,n);A=rand(n,n);
+x=rdiv(b,A);
+if Err(x*A-b)>10000*%eps then bugmes();quit;end
+//Small dimensions complex
+b=rand(2,n)+%i;A=rand(n,n);
+x=rdiv(b,A);
+if Err(x*A-b)>10000*%eps then bugmes();quit;end
+b=rand(2,n);A=rand(n,n)+%i;
+x=rdiv(b,A);
+if Err(x*A-b)>10000*%eps then bugmes();quit;end
+b=rand(2,n)+%i;A=rand(n,n)+%i;
+x=rdiv(b,A);
+if Err(x*A-b)>10000*%eps then bugmes();quit;end
+//-----Rectangular matrix-----
+n=5;m=3; A=rand(m,n);b=rand(2,n+1);
+if rdiv([],A) <>[] then bugmes();quit;end
+if execstr('rdiv(b,A)','errcatch')==0 then bugmes();quit;end
+//Small dimensions real
+n=5;m=3;
+b=rand(2,n);A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>200*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(2,n);A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>200*%eps then bugmes();quit;end
+//Small dimensions complex
+n=5;m=3;
+b=rand(2,n)+%i;A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>200*%eps then bugmes();quit;end
+n=5;m=3;
+b=rand(2,n);A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>200*%eps then bugmes();quit;end
+b=rand(2,n)+%i;A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>200*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(2,n)+%i;A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(2,n);A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=3;m=5;
+b=rand(2,n)+%i;A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+//LArge dimension real
+n=50;m=30;
+b=rand(2,n);A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=30;m=50;
+b=rand(2,n);A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+//Large dimensions complex
+n=50;m=30;
+b=rand(2,n)+%i;A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=50;m=30;
+b=rand(2,n);A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+b=rand(2,n)+%i;A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=30;m=50;
+b=rand(2,n)+%i;A=rand(m,n);
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=30;m=50;
+b=rand(2,n);A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
+n=30;m=50;
+b=rand(2,n)+%i;A=rand(m,n)+%i;
+x=rdiv(b,A);
+if Err(x*A*A'-b*A')>1000*%eps then bugmes();quit;end
index 504bfce..6e558dd 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')
 
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/schur.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/schur.dia.ref
new file mode 100644 (file)
index 0000000..50a543a
--- /dev/null
@@ -0,0 +1,407 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x)
+       r=norm(x,1)
+endfunction
+rand('normal')
+//define tools
+function A=testmat1(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+       A(1,n)=1/a;A(n,1)=a
+endfunction
+//==========================================================================
+//==============================    schur     ==============================
+//==========================================================================
+clear sel
+function t=sel(R),t=real(R)<0 ,endfunction
+//Empty matrix
+A=[];
+if schur(A)<>[] then bugmes();quit;end
+if schur(A,'real')<>[] then bugmes();quit;end
+if schur(A,'complex')<>[] then bugmes();quit;end
+if schur(A,'c')<>[] then bugmes();quit;end
+if schur(A,'d')<>[] then bugmes();quit;end
+if schur(A,sel)<>[] then bugmes();quit;end
+[U,S]=schur(A);
+if U<>[]|S<>[] then bugmes();quit;end
+[U,S]=schur(A,'real');
+if U<>[]|S<>[] then bugmes();quit;end
+[U,S]=schur(A,'complex');
+if U<>[]|S<>[] then bugmes();quit;end
+[U,N]=schur(A,'c');
+if U<>[]|N<>0 then bugmes();quit;end
+[U,N]=schur(A,'d');
+if U<>[]|N<>0 then bugmes();quit;end
+[U,N]=schur(A,sel);
+if U<>[]|N<>0 then bugmes();quit;end
+[U,N,S]=schur(A,'c');
+if U<>[]|N<>0|S<>[] then bugmes();quit;end
+[U,N,S]=schur(A,'d');
+if U<>[]|N<>0|S<>[] then bugmes();quit;end
+[U,N,S]=schur(A,sel);
+if U<>[]|N<>0|S<>[] then bugmes();quit;end
+//Rectangular matrix
+if execstr('schur(rand(2,3))','errcatch')==0 then bugmes();quit;end
+if execstr('[U,S]=schur(rand(2,3))','errcatch')==0 then bugmes();quit;end
+if execstr('schur(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
+if execstr('[U,S]=schur(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
+//Small dimension
+A=testmat1(3,5);Ac=testmat1(3+%i,5);
+//Real
+[U,S]=schur(A);
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
+if Err(schur(A)-S) >%eps then bugmes();quit;end
+[U,S]=schur(A,'real');
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
+if Err(schur(A)-S) >%eps then bugmes();quit;end
+[U,S]=schur(A,'complex');
+if Err(triu(S)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
+if Err(schur(A,'complex')-S) >%eps then bugmes();quit;end
+[U,n]=schur(A,'c');S=U'*A*U;
+if n<>2 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+[U,n]=schur(A,'d');S=U'*A*U;
+if n<>0 then bugmes();quit;end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
+[U,n]=schur(A,sel);S=U'*A*U;
+if n<>2 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+//Complex
+[U,S]=schur(Ac);
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-Ac)>200*%eps then bugmes();quit;end
+if Err(schur(Ac)-S) >%eps then bugmes();quit;end
+[U,S]=schur(Ac,'complex');
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-Ac)>200*%eps then bugmes();quit;end
+if Err(schur(Ac)-S) >%eps then bugmes();quit;end
+[U,n]=schur(Ac,'c');S=U'*Ac*U;
+if n<>3 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+[U,n]=schur(Ac,'d');S=U'*A*U;
+if n<>0 then bugmes();quit;end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
+[U,n]=schur(Ac,sel);S=U'*Ac*U;
+if n<>3 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+//Large dimension
+A=testmat1(3,50);Ac=testmat1(3+%i,50);
+//Real
+[U,S]=schur(A);
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
+if Err(schur(A)-S) >%eps then bugmes();quit;end
+[U,S]=schur(A,'real');
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
+if Err(schur(A)-S) >%eps then bugmes();quit;end
+[U,S]=schur(A,'complex');
+if Err(triu(S)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
+if Err(schur(A,'complex')-S) >%eps then bugmes();quit;end
+[U,n]=schur(A,'c');S=U'*A*U;
+if n<>25 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+[U,n]=schur(A,'d');S=U'*A*U;
+if n<>0 then bugmes();quit;end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
+[U,n]=schur(A,sel);S=U'*A*U;
+if n<>25 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+//Complex
+[U,S]=schur(Ac);
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-Ac)>1000*%eps then bugmes();quit;end
+if Err(schur(Ac)-S) >%eps then bugmes();quit;end
+[U,S]=schur(Ac,'complex');
+if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
+if Err(U*S*U'-Ac)>1000*%eps then bugmes();quit;end
+if Err(schur(Ac)-S) >%eps then bugmes();quit;end
+[U,n]=schur(Ac,'c');S=U'*Ac*U;
+if n<>25 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+[U,n]=schur(Ac,'d');S=U'*Ac*U;
+if n<>0 then bugmes();quit;end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
+[U,n]=schur(Ac,sel);S=U'*Ac*U;
+if n<>25 then bugmes();quit;end
+if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
+//==========================================================================
+//==============================    schur part II   ========================
+//==========================================================================
+//Empty matrix
+[As,Es]=schur([],[]);
+if As<>[]|Es<>[] then bugmes();quit;end
+[As,dim]=schur([],[],'c');
+if As<>[]|dim<>0 then bugmes();quit;end
+[As,dim]=schur([],[],'d');
+if As<>[]|dim<>0 then bugmes();quit;end
+[As,dim]=schur([],[],sel);
+if As<>[]|dim<>0 then bugmes();quit;end
+[As,Es,Q,Z]=schur([],[]);
+if As<>[]|Es<>[]|Q<>[]|Z<>[] then bugmes();quit;end
+[As,Es,dim]=schur([],[],'c');
+if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
+[As,Es,dim]=schur([],[],'d');
+if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
+[As,Es,dim]=schur([],[],sel);
+if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
+[Z,dim]=schur([],[],'c');
+if Z<>[]|dim<>0 then bugmes();quit;end
+[Z,dim]=schur([],[],'d');
+if Z<>[]|dim<>0 then bugmes();quit;end
+[Z,dim]=schur([],[],sel);
+if Z<>[]|dim<>0 then bugmes();quit;end
+//Rectangular matrix
+if execstr('[As,Es]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  bugmes();quit;end
+if execstr('[As,Es,Q,Z]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  bugmes();quit;end
+if execstr('[As,Es,dim]=schur(rand(2,3),rand(2,3),''c'')','errcatch')==0 then  bugmes();quit;end
+if execstr('[Z,dim]=schur(rand(2,3),rand(2,3),sel)','errcatch')==0 then  bugmes();quit;end
+//Small dimension
+//----Real------------
+A=testmat1(1,5);E=testmat1(-2,5) ;
+[As,Es,Q,Z]=schur(A,E);
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+[As1,Es1]=schur(A,E);
+if Err(As1-As)>10*%eps then bugmes();quit;end
+if Err(Es1-Es)>10*%eps then bugmes();quit;end
+// Ordered 'c'
+dim=schur(A,E,'c');
+if dim<>5 then bugmes();quit;end
+[Z,dim]=schur(A,E,'c');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'c');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>5 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+// Ordered 'd'
+dim=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+[Z,dim]=schur(A,E,'d');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'d');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>5 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//ordered sel
+clear sel
+function t=sel(Alpha,Beta),t=real(Alpha)>-0.2*real(Beta) ,endfunction
+dim=schur(A,E,sel);
+if dim<>2 then bugmes();quit;end
+[Z,dim]=schur(A,E,sel);
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,sel);
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>2 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,sel);
+if dim<>2 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//----Complex------------
+A=testmat1(1+%i,5);E=testmat1(-2-3*%i,5) ;
+[As,Es,Q,Z]=schur(A,E);
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+[As1,Es1]=schur(A,E);
+if Err(As1-As)>10*%eps then bugmes();quit;end
+if Err(Es1-Es)>10*%eps then bugmes();quit;end
+// Ordered 'c'
+dim=schur(A,E,'c');
+if dim<>5 then bugmes();quit;end
+[Z,dim]=schur(A,E,'c');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'c');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>5 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+// Ordered 'd'
+dim=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+[Z,dim]=schur(A,E,'d');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'d');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>5 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>5 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//ordered sel
+clear sel
+function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
+dim=schur(A,E,sel);
+if dim<>3 then bugmes();quit;end
+[Z,dim]=schur(A,E,sel);
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,sel);
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>3 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,sel);
+if dim<>3 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//Large dimension
+//----Real------------
+A=testmat1(1,50);E=testmat1(-2,50) ;
+[As,Es,Q,Z]=schur(A,E);
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+[As1,Es1]=schur(A,E);
+if Err(As1-As)>10*%eps then bugmes();quit;end
+if Err(Es1-Es)>10*%eps then bugmes();quit;end
+// Ordered 'c'
+dim=schur(A,E,'c');
+if dim<>50 then bugmes();quit;end
+[Z,dim]=schur(A,E,'c');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'c');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>50 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+// Ordered 'd'
+dim=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+[Z,dim]=schur(A,E,'d');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'d');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>50 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//ordered sel
+clear sel
+function t=sel(Alpha,Beta)
+       t=real(Alpha)>-0.2*real(Beta)
+endfunction
+dim=schur(A,E,sel); // plante ici DGGES LAPACK 3.1
+if dim<>12 then bugmes();quit;end
+[Z,dim]=schur(A,E,sel);
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,sel);
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>12 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,sel);
+if dim<>12 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
+//----Complex------------
+A=testmat1(1+%i,50);E=testmat1(-2-3*%i,50) ;
+[As,Es,Q,Z]=schur(A,E);
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
+[As1,Es1]=schur(A,E);
+if Err(As1-As)>10*%eps then bugmes();quit;end
+if Err(Es1-Es)>10*%eps then bugmes();quit;end
+// Ordered 'c'
+dim=schur(A,E,'c');
+if dim<>50 then bugmes();quit;end
+[Z,dim]=schur(A,E,'c');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'c');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>50 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
+// Ordered 'd'
+dim=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+[Z,dim]=schur(A,E,'d');
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,'d');
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>50 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,'d');
+if dim<>50 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
+//ordered sel
+clear sel
+function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
+dim=schur(A,E,sel);
+if dim<>32 then bugmes();quit;end
+[Z,dim]=schur(A,E,sel);
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+[Q,Z1,dim]=schur(A,E,sel);
+if Err(Z1-Z)>10*%eps then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if dim<>32 then bugmes();quit;end
+[As,Es,Z,dim]=schur(A,E,sel);
+if dim<>32 then bugmes();quit;end
+if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
+if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
+if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
+if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
index f6a1607..1e3fc34 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x)
        r=norm(x,1)
 endfunction
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/spec.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/spec.dia.ref
new file mode 100644 (file)
index 0000000..2328781
--- /dev/null
@@ -0,0 +1,86 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//define tools
+function A=testmat1(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
+       A(1,n)=1/a;A(n,1)=a
+endfunction
+function r=Err(x)
+       r=norm(x,1)
+endfunction
+rand('normal')
+//==========================================================================
+//==============================     spec     ==============================
+//==========================================================================
+function r=Checktestmat1(a,n)
+   A=testmat1(a,n);S=spec(A);
+   SR=real(S);SI=imag(S);
+   dt=2*%i*%pi/n;Z=exp(dt*(1:n)');ZR=real(Z*((1+a*a')/a));
+   ZI=-imag(Z*((a*a'-1)/a));
+   r=max(norm(sort(SR)-sort(ZR)),norm(sort(SI)-sort(ZI)))
+endfunction
+function A=testmat2(a,n)
+       //eigen values are given by a dilation of nth roots of 1
+       A=testmat1(a,n);A=A+A'
+endfunction
+function r=Checktestmat2(a,n)
+   A=testmat2(a,n);S=spec(A);
+   SR=real(S);SI=imag(S);
+   dt=2*%i*%pi/n;Z=exp(dt*(1:n)');
+   ZR=2*real(Z*((1+a*a')/a));ZI=0*ZR;
+   r=max(norm(sort(SR)-sort(ZR)),norm(sort(SI)-sort(ZI)))
+endfunction
+//Empty matrix
+A=[];
+S=spec(A);
+if S<>[] then bugmes();quit;end
+//Matrix with Inf or Nan (test de la detection d'erreur
+if execstr('spec([%inf 1;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('spec([1 %nan;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('spec([%inf %i;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('spec([%i %i;%nan 3])','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+//Real Case
+//Unsymetric
+if Checktestmat1(3,5)>200*%eps then bugmes();quit;end
+[U,S]=spec(testmat1(3,5));
+if Err(U*S/U-testmat1(3,5))>200*%eps then bugmes();quit;end
+//Symmetric
+if Checktestmat2(3,5)>200*%eps then bugmes();quit;end
+[U,S]=spec(testmat2(3,5));
+if Err(U*S/U-testmat2(3,5))>200*%eps then bugmes();quit;end
+//Complex Case
+//Unsymetric
+if Checktestmat1(3+2*%i,5)>200*%eps then bugmes();quit;end
+[U,S]=spec(testmat1(3+2*%i,5));
+if Err(U*S/U-testmat1(3+2*%i,5))>200*%eps then bugmes();quit;end
+//Symmetric
+if Checktestmat2(3+2*%i,5)>200*%eps then bugmes();quit;end
+[U,S]=spec(testmat2(3+2*%i,5));
+if Err(U*S/U-testmat2(3+2*%i,5))>200*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+//Real Case
+//Unsymetric
+if Checktestmat1(3,50)>1000*%eps then bugmes();quit;end
+[U,S]=spec(testmat1(3,50));
+if Err(U*S/U-testmat1(3,50))>1000*%eps then bugmes();quit;end
+//Symmetric
+if Checktestmat2(3,50)>1000*%eps then bugmes();quit;end
+[U,S]=spec(testmat2(3,50));
+if Err(U*S/U-testmat2(3,50))>1000*%eps then bugmes();quit;end
+//Complex Case
+//Unsymetric
+if Checktestmat1(3+2*%i,50)>1000*%eps then bugmes();quit;end
+[U,S]=spec(testmat1(3+2*%i,50));
+if Err(U*S/U-testmat1(3+2*%i,50))>1000*%eps then bugmes();quit;end
+//Symmetric
+if Checktestmat2(3+2*%i,50)>1000*%eps then bugmes();quit;end
+[U,S]=spec(testmat2(3+2*%i,50));
+if Err(U*S/U-testmat2(3+2*%i,50))>1000*%eps then bugmes();quit;end
index 01928c4..fd50506 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 //define tools
 function A=testmat1(a,n)
        //eigen values are given by a dilation of nth roots of 1
diff --git a/scilab/modules/linear_algebra/tests/unit_tests/svd.dia.ref b/scilab/modules/linear_algebra/tests/unit_tests/svd.dia.ref
new file mode 100644 (file)
index 0000000..f324882
--- /dev/null
@@ -0,0 +1,300 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+function r=Err(x),r=norm(x,1),endfunction
+rand('normal')
+//==========================================================================
+//==============================     svd      ==============================
+//==========================================================================
+//Empty matrix
+A=[];
+[U,S,V]=svd(A);
+if U<>[]|V<>[]|S<>[] then bugmes();quit;end
+S=svd(A);
+if S<>[] then bugmes();quit;end
+[U,S,V]=svd(A,"e");
+if U<>[]|V<>[]|S<>[] then bugmes();quit;end
+S=svd(A,"e");
+if S<>[] then bugmes();quit;end
+//Matrix with inf or nan
+if execstr('svd([%inf 1;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('svd([1 %nan;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('svd([%inf %i;2 3])','errcatch')==0 then bugmes();quit;end
+if execstr('svd([%i %i;%nan 3])','errcatch')==0 then bugmes();quit;end
+//Small dimension
+//---------------
+A=rand(3,5);Ac=A+%i*rand(A);
+//Real Case
+[U,S,V]=svd(A);
+if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
+if Err(svd(A)-diag(S))> 200*%eps then bugmes();quit;end
+[U,S,V]=svd(A,"e");
+if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
+A=A';
+[U,S,V]=svd(A);
+if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
+if Err(svd(A)-diag(S))> 200*%eps then bugmes();quit;end
+[U,S,V]=svd(A,"e");
+if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
+//Complex Case
+[U,S,V]=svd(Ac);
+if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
+if Err(svd(Ac)-diag(S))> 200*%eps then bugmes();quit;end
+[U,S,V]=svd(Ac,"e");
+if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
+Ac=Ac';
+[U,S,V]=svd(Ac);U*S*V'-Ac;
+if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
+if Err(svd(Ac)-diag(S))> 200*%eps then bugmes();quit;end
+[U,S,V]=svd(Ac,"e");
+if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
+//Large dimension
+//---------------
+A=rand(150,60);Ac=A+rand(A)*%i;
+//Real Case
+[U,S,V]=svd(A);
+if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
+if Err(svd(A)-diag(S))> 10000*%eps then bugmes();quit;end
+[U,S,V]=svd(A,"e");
+if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
+A=A';
+[U,S,V]=svd(A);
+if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
+if Err(svd(A)-diag(S))> 10000*%eps then bugmes();quit;end
+[U,S,V]=svd(A,"e");
+if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
+//Complex Case
+[U,S,V]=svd(Ac);
+if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
+if Err(svd(Ac)-diag(S))> 10000*%eps then bugmes();quit;end
+[U,S,V]=svd(Ac,"e");
+if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
+Ac=Ac';
+[U,S,V]=svd(Ac);U*S*V'-Ac;
+if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
+if Err(svd(Ac)-diag(S))> 10000*%eps then bugmes();quit;end
+[U,S,V]=svd(Ac,"e");
+if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
+//==========================================================================
+//==============================     svd part II     =======================
+//==========================================================================
+//Empty matrix
+if svd([])<>[] then bugmes();quit;end
+if svd([],"e")<>[] then bugmes();quit;end
+[U,S]=svd([])
+ S  =
+     []
+ U  =
+     []
+if U<>[]|S<>[]  then bugmes();quit;end
+[U,S,V]=svd([]);
+if U<>[]|S<>[]|V<>[]  then bugmes();quit;end
+[U,S,V,rk]=svd([]);
+if U<>[]|S<>[]|V<>[]|rk<>0  then bugmes();quit;end
+[U,S,V,rk]=svd([],%eps);
+if U<>[]|S<>[]|V<>[]|rk<>0  then bugmes();quit;end
+if execstr('[U,S,V,rk]=svd([],'"e'")','errcatch') == 0 then bugmes();quit;end
+//Small dimension
+//Real Case Fat
+A=rand(3,5);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>10*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>10*%eps  then bugmes();quit;end
+if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>200*%eps  then bugmes();quit;end
+if Err(V-V1) >200*%eps  then bugmes();quit;end
+if rk<>3 then bugmes();quit;end
+//Real Case Tall
+A=rand(5,3);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if size(U1,2)<>3 then bugmes();quit;end
+if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>200*%eps  then bugmes();quit;end
+if Err(V-V1) >200*%eps  then bugmes();quit;end
+if rk<>3 then bugmes();quit;end
+//Complex Case Fat
+A=rand(3,5)+%i*rand(3,5);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >30*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>200*%eps  then bugmes();quit;end
+if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>200*%eps  then bugmes();quit;end
+if Err(V-V1) >200*%eps  then bugmes();quit;end
+if rk<>3 then bugmes();quit;end
+//Complex Case Tall
+A=rand(5,3)+%i*rand(5,3);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if size(U1,2)<>3 then bugmes();quit;end
+if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
+if Err(U-U1)>200*%eps  then bugmes();quit;end
+if Err(V-V1) >200*%eps  then bugmes();quit;end
+if rk<>3 then bugmes();quit;end
+//Large dimension
+//Real Case Fat
+A=rand(30,50);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>1000*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>1000*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>1000*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >1000*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
+if Err(U-U1)>10*%eps  then bugmes();quit;end
+if Err(U1*S1*V1'-A) >1000*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
+if Err(U-U1)>5000*%eps  then bugmes();quit;end
+if Err(V-V1) >5000*%eps  then bugmes();quit;end
+if rk<>30 then bugmes();quit;end
+//Real Case Tall
+A=rand(50,30);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if size(U1,2)<>30 then bugmes();quit;end
+if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U-U1)>5000*%eps  then bugmes();quit;end
+if Err(V-V1) >5000*%eps  then bugmes();quit;end
+if rk<>30 then bugmes();quit;end
+//Complex Case Fat
+A=rand(30,50)+%i*rand(30,50);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U-U1)>5000*%eps  then bugmes();quit;end
+if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U-U1)>5000*%eps  then bugmes();quit;end
+if Err(V-V1) >5000*%eps  then bugmes();quit;end
+if rk<>30 then bugmes();quit;end
+//Complex Case Tall
+A=rand(50,30)+%i*rand(50,30);
+S=svd(A);
+if or(S<0) then bugmes();quit;end
+if sort(S)<>S  then bugmes();quit;end
+[U,S1]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
+[U1,S1,V]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
+if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1]=svd(A,"e");
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if size(U1,2)<>30 then bugmes();quit;end
+if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
+[U1,S1,V1,rk]=svd(A);
+if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
+if Err(U-U1)>5000*%eps  then bugmes();quit;end
+if Err(V-V1) >5000*%eps  then bugmes();quit;end
+if rk<>30 then bugmes();quit;end
+function c=cond(A)
+  if A==[] then c=1,else S=svd(A);c=S($)/S(1),end
+endfunction
+Warning : redefining function: cond                    . Use funcprot(0) to avoid this message
+
index 5381667..b6712a2 100644 (file)
@@ -4,7 +4,6 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-// <-- NO CHECK REF -->
 function r=Err(x),r=norm(x,1),endfunction
 rand('normal')