linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / schur.tst
index 1e3fc34..bd6c3ae 100644 (file)
@@ -8,6 +8,12 @@ function r=Err(x)
        r=norm(x,1)
 endfunction
 rand('normal')
+ilib_verbose(0);
+
+libName = [];
+if getos() == "Windows" then
+    libName = SCI + "/bin/elem_func";
+end
 
 //define tools
 function A=testmat1(a,n)
@@ -178,6 +184,83 @@ if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
 if n<>25 then pause,end
 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
+
+// Lib part
+cd TMPDIR;
+// equal to schur(A, 'c');
+C=[ 'int mytest1(double* _real, double* _img)'
+    '{' 
+    '    return *_real < 0;'
+    '}'];
+
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest1','mytest.c',[],'c');
+exec loader.sce;
+[U,n]=schur(A,'mytest1');S=U'*A*U;
+if n<>25 then pause,end
+if or(real(spec(S(1:n,1:n)))>=0) then pause,end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
+
+// equal to schur(A, 'd');
+C=[ 'extern double dpythags(double,double);' // YaSp function
+    ''
+    'int mytest2(double* _real, double* _img)'
+    '{' 
+    '    return dpythags(*_real, *_img) < 1;'
+    '}'];
+
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest2','mytest.c', libName,'c');
+
+exec loader.sce;
+[U,n]=schur(A,'mytest2');S=U'*A*U;
+if n<>0 then pause,end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
+
+// equal to schur(Ac, 'c');
+C=[ '#include ""doublecomplex.h""'
+    ''
+    'int mytest3(doublecomplex* _complex)'
+    '{'
+    '    return _complex->r < 0 ? 1 : 0;'
+    '}'];
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest3','mytest.c',[],'c');
+exec loader.sce;
+[U,n]=schur(Ac, 'mytest3');S=U'*Ac*U;
+if n<>25 then pause,end
+if or(real(spec(S(1:n,1:n)))>=0) then pause,end
+if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
+
+// equal to schur(Ac, 'd');
+C=[ '#include ""doublecomplex.h""'
+    ''
+    'extern double dpythags(double,double);' // YaSp function
+    ''
+    'int mytest4(doublecomplex* _complex)'
+    '{'
+    '    if(dpythags(_complex->r, _complex->i) < 1)'
+    '    {'
+    '        return 1;'
+    '    }'
+    '    else'
+    '    {'
+    '        return 0;'
+    '    }'
+    '}'];
+
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest4','mytest.c', libName,'c');
+exec loader.sce;
+[U,n]=schur(Ac,'mytest4');S=U'*Ac*U;
+if n<>0 then pause,end
+if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
+
+
 //==========================================================================
 //==============================    schur part II   ======================== 
 //==========================================================================
@@ -240,7 +323,7 @@ if Err(Z1-Z)>10*%eps then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if dim<>5 then pause,end
 
-[As,Es,Z,dim]=schur(A,E,'d');
+[As,Es,Z,dim]=schur(A,E,'c');
 if dim<>5 then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
@@ -266,7 +349,9 @@ if Err(Es-Q'*E*Z) >200*%eps then pause,end
 
 //ordered sel
 clear sel
-function t=sel(Alpha,Beta),t=real(Alpha)>-0.2*real(Beta) ,endfunction
+function t=sel(Alpha,Beta)
+    t = real(Alpha) > -0.2 * real(Beta);
+endfunction
 
 dim=schur(A,E,sel);
 if dim<>2 then pause,end
@@ -284,6 +369,8 @@ if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
 if Err(As-Q'*A*Z) >200*%eps then pause,end
 if Err(Es-Q'*E*Z) >200*%eps then pause,end
+
+
 //----Complex------------
 A=testmat1(1+%i,5);E=testmat1(-2-3*%i,5) ;
 [As,Es,Q,Z]=schur(A,E);
@@ -307,12 +394,14 @@ if Err(Z1-Z)>10*%eps then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if dim<>5 then pause,end
 
-[As,Es,Z,dim]=schur(A,E,'d');
+[As,Es,Z,dim]=schur(A,E,'c');
 if dim<>5 then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
 if Err(As-Q'*A*Z) >200*%eps then pause,end
 if Err(Es-Q'*E*Z) >200*%eps then pause,end
+
+
 // Ordered 'd'
 dim=schur(A,E,'d');
 if dim<>5 then pause,end
@@ -377,7 +466,7 @@ if Err(Z1-Z)>10*%eps then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if dim<>50 then pause,end
 
-[As,Es,Z,dim]=schur(A,E,'d');
+[As,Es,Z,dim]=schur(A,E,'c');
 if dim<>50 then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
@@ -404,10 +493,10 @@ if Err(Es-Q'*E*Z) >200*%eps then pause,end
 //ordered sel
 clear sel
 function t=sel(Alpha,Beta)
-       t=real(Alpha)>-0.2*real(Beta)
+       t = real(Alpha) > -0.2 * real(Beta);
 endfunction
 
-dim=schur(A,E,sel); // plante ici DGGES LAPACK 3.1
+dim=schur(A,E,sel);
 if dim<>12 then pause,end
 [Z,dim]=schur(A,E,sel);
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
@@ -423,6 +512,76 @@ if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
 if Err(As-Q'*A*Z) >200*%eps then pause,end
 if Err(Es-Q'*E*Z) >200*%eps then pause,end
+
+// Lib part
+cd TMPDIR;
+
+// equal to schur(A, E, 'c');
+C=[ 'extern double dlamch_(char *CMACH, unsigned long int);' // dlamch_ == C2F(dlamch)
+    ''
+    'int mytest5(double* _real, double* _img, double* _beta)'
+    '{'
+    '    double dblP = dlamch_((char*)""p"", 1L);'
+    '    int iTest1 = (*_real < 0 && *_beta > 0);'
+    '    int iTest2 = (*_real > 0 && *_beta < 0);'
+    '    int iTest3 = (fabs(*_beta) > fabs(*_real) * dblP);'
+    ''
+    '    return (iTest1 || iTest2 && iTest3);'
+    '}'];
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest5','mytest.c',[],'c');
+exec loader.sce;
+
+dim=schur(A,E,'mytest5');
+if dim<>50 then pause,end
+[Z,dim]=schur(A,E,'mytest5');
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+
+[Q,Z1,dim]=schur(A,E,'mytest5');
+if Err(Z1-Z)>10*%eps then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if dim<>50 then pause,end
+
+[As,Es,Z,dim]=schur(A,E,'mytest5');
+if dim<>50 then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+if Err(As-Q'*A*Z) >200*%eps then pause,end
+if Err(Es-Q'*E*Z) >200*%eps then pause,end
+
+// equal to schur(A, E, 'd');
+C=[ 'extern double dpythags(double,double);' // YaSp function
+    ''
+    'int mytest6(double* _real, double* _img, double* _beta)'
+    '{'
+    '    double dblPythag =  dpythags(*_real, *_img);'
+    ''
+    '    return (dblPythag < fabs(*_beta));'
+    '}'];
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest6','mytest.c', libName,'c');
+exec loader.sce;
+
+dim=schur(A,E,'mytest6');
+if dim<>50 then pause,end
+
+[Z,dim]=schur(A,E,'mytest6');
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+
+[Q,Z1,dim]=schur(A,E,'mytest6');
+if Err(Z1-Z)>10*%eps then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if dim<>50 then pause,end
+
+[As,Es,Z,dim]=schur(A,E,'mytest6');
+if dim<>50 then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+if Err(As-Q'*A*Z) >200*%eps then pause,end
+if Err(Es-Q'*E*Z) >200*%eps then pause,end
+
 //----Complex------------
 A=testmat1(1+%i,50);E=testmat1(-2-3*%i,50) ;
 [As,Es,Q,Z]=schur(A,E);
@@ -446,7 +605,7 @@ if Err(Z1-Z)>10*%eps then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if dim<>50 then pause,end
 
-[As,Es,Z,dim]=schur(A,E,'d');
+[As,Es,Z,dim]=schur(A,E,'c');
 if dim<>50 then pause,end
 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
@@ -491,3 +650,69 @@ if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
 if Err(As-Q'*A*Z) >1000*%eps then pause,end
 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
 
+//Lib part
+// equal to schur(A, E, 'c');
+C=[ '#include ""doublecomplex.h"";'
+    ''
+    'int mytest7(doublecomplex* _complex)'
+    '{'
+    '    return (_complex->r < 0);'
+    '}'];
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest7','mytest.c',[],'c');
+exec loader.sce;
+
+dim=schur(A,E,'mytest7');
+if dim<>50 then pause,end
+
+[Z,dim]=schur(A,E,'mytest7');
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+
+[Q,Z1,dim]=schur(A,E,'mytest7');
+if Err(Z1-Z)>10*%eps then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if dim<>50 then pause,end
+
+[As,Es,Z,dim]=schur(A,E,'mytest7');
+if dim<>50 then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+if Err(As-Q'*A*Z) >1000*%eps then pause,end
+if Err(Es-Q'*E*Z) >1000*%eps then pause,end
+
+// equal to schur(A, E, 'd');
+C=[ '#include ""doublecomplex.h"";'
+    ''
+    'extern double dpythags(double,double);' // YaSp function
+    ''
+    'int mytest8(doublecomplex* _alpha, doublecomplex* _beta)'
+    '{'
+    '    double dblP1 = dpythags(_alpha->r, _alpha->i);'
+    '    double dblP2 = dpythags(_beta->r, _beta->i);'
+    ''
+    '    return (dblP1 <  dblP2);'
+    '}'];
+mputl(C,TMPDIR+'/mytest.c');
+ulink();
+lp=ilib_for_link('mytest8','mytest.c', libName,'c');
+exec loader.sce;
+
+dim=schur(A,E,'mytest8');
+if dim<>50 then pause,end
+
+[Z,dim]=schur(A,E,'mytest8');
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+
+[Q,Z1,dim]=schur(A,E,'mytest8');
+if Err(Z1-Z)>10*%eps then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if dim<>50 then pause,end
+
+[As,Es,Z,dim]=schur(A,E,'mytest8');
+if dim<>50 then pause,end
+if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
+if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
+if Err(As-Q'*A*Z) >1000*%eps then pause,end
+if Err(Es-Q'*E*Z) >1000*%eps then pause,end
+