* Bug 15701 fixed: now A\B is faster when A is square and triangular 70/20970/6
Stéphane MOTTELET [Tue, 7 May 2019 15:44:14 +0000 (17:44 +0200)]
http://bugzilla.scilab.org/show_bug.cgi?id=15701

Change-Id: Ib311ac1cdd0dbc98439daa23c50231d66c29da6e

scilab/CHANGES.md
scilab/modules/ast/includes/operations/matrix_division.h
scilab/modules/ast/src/c/operations/matrix_division.c
scilab/modules/ast/tests/nonreg_tests/bug_15701.tst [new file with mode: 0644]
scilab/modules/elementary_functions/includes/elem_common.h

index d2d3595..0612862 100644 (file)
@@ -255,6 +255,7 @@ Bug Fixes
 * [#15523](http://bugzilla.scilab.org/show_bug.cgi?id=15523): `%ODEOPTIONS(1)=2` didn't work with solvers 'rk' and 'rkf'
 * [#15534](http://bugzilla.scilab.org/show_bug.cgi?id=15534): Booleans and encoded integers could not be concatenated together.
 * [#15577](http://bugzilla.scilab.org/show_bug.cgi?id=15577): `edit` did not accept a line number as text, as with `edit linspace 21`.
+
 * [#15580](http://bugzilla.scilab.org/show_bug.cgi?id=15580): `det(sparse([],[]))` yielded an error.
 * [#15981](http://bugzilla.scilab.org/show_bug.cgi?id=15981): `wavread()` kept the wav file open and locked when returning on errors. It weakly managed the input file name. It claimed for invalid data formats instead of unsupported ones, with poor information about the current format vs the supported ones. Several error messages refered to a wrong function.
 * [#15595](http://bugzilla.scilab.org/show_bug.cgi?id=15595): `unique()` was not able to return distinct values in their original order, without sorting them. A `keepOrder` option now allows it.
@@ -263,6 +264,7 @@ Bug Fixes
 * [#15742](http://bugzilla.scilab.org/show_bug.cgi?id=15742): The `compatibility_functions` module should be merged in the `m2sci` one.
 * [#15581](http://bugzilla.scilab.org/show_bug.cgi?id=15581): display of complex matrix was ugly.
 * [#15680](http://bugzilla.scilab.org/show_bug.cgi?id=15680): `loadmatfile` could not return variables in a structure instead of into the calling environment.
+* [#15701](http://bugzilla.scilab.org/show_bug.cgi?id=15701): `A\B` was not faster when `A` is square and triangular.
 * [#15734](http://bugzilla.scilab.org/show_bug.cgi?id=15734):  Trivial infinite loop could not be interrupted.
 * [#15744](http://bugzilla.scilab.org/show_bug.cgi?id=15744): `sylm(a,b)` yielded an error when degree(a)==0 or degree(b)==0.
 * [#15745](http://bugzilla.scilab.org/show_bug.cgi?id=15745): `diophant(0,0,m)`, `diophant([p 0],q)`, `diophant([0 p],q)` with m<>0 and p>q were wrong. There was no flag for cases with an infinite number of solutions. When there is no solution, some values were returned anyway, instead of []. In this case, the documented definition of the err value was dubious. Decimal numbers and integers were accepted, but not encoded integers. Inf and NaN input coefficients were not rejected.
index 3189fbb..01bea5e 100644 (file)
@@ -22,4 +22,7 @@
 #include "matrix_left_division.h"
 #include "matrix_right_division.h"
 
+void matrixIsTriangular(double *_pdblReal, double *_pdblImg, int _iRows, int _iCols, int *bUpper, int *bLower);
+extern double C2F(dasum) (int *, double *, int *);
+
 #endif //!__MATRIX_DIV_H__
index 70e4dd1..9c70d9e 100644 (file)
@@ -1,6 +1,7 @@
 /*
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008-2008 - INRIA - Antoine ELIAS <antoine.elias@scilab.org>
+ * Copyright (C) 2019 - Stéphane MOTTELET
  *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
@@ -565,6 +566,8 @@ int iLeftDivisionOfRealMatrix(
     int iIndex  = 0;
     char cNorm  = 0;
     int iExit   = 0;
+    int iLower = 0;
+    int iUpper = 0;
 
     /*temporary variables*/
     int iWorkMin    = 0;
@@ -575,6 +578,7 @@ int iLeftDivisionOfRealMatrix(
     double dblEps       = 0;
     double RCONDthresh  = 0;
     double dblAnorm     = 0;
+    double dblOne       = 1;
 
     double *pAf     = NULL;
     double *pXb     = NULL;
@@ -600,37 +604,58 @@ int iLeftDivisionOfRealMatrix(
     pJpvt       = (int*)malloc(sizeof(int) * _iCols1);
     pIwork      = (int*)malloc(sizeof(int) * _iCols1);
 
-    cNorm       = '1';
     pDwork      = (double*)malloc(sizeof(double) * iWorkMin);
     dblEps      = nc_eps();
     RCONDthresh = 10 * dblEps;
-    dblAnorm    = C2F(dlange)(&cNorm, &_iRows1, &_iCols1, _pdblReal1, &_iRows1, pDwork);
 
     if (_iRows1 == _iCols1)
     {
-        cNorm = 'F';
-        C2F(dlacpy)(&cNorm, &_iCols1, &_iCols1, _pdblReal1, &_iCols1, pAf, &_iCols1);
-        C2F(dgetrf)(&_iCols1, &_iCols1, pAf, &_iCols1, pIpiv, &iInfo);
-        if (iInfo == 0)
+        matrixIsTriangular(_pdblReal1, NULL, _iRows1, _iCols1, &iUpper, &iLower);
+        if (iUpper || iLower)
         {
-            cNorm = '1';
-            C2F(dgecon)(&cNorm, &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
+            //if matrix is triangular
+            char cUpLoType[4] = {'N','U','L','U'};
+            char cUpLo = cUpLoType[iUpper + 2*iLower];
+            C2F(dtrcon)("1", &cUpLo, "N", &_iRows1, _pdblReal1, &_iRows1, &dblRcond, pDwork, pIwork, &iInfo);
             if (dblRcond > RCONDthresh)
             {
-                // _pdblReal2 will be overwrite by dgetrs
+                // _pdblReal2 will be overwritten by dtrsm
                 iSize = _iRows2 * _iCols2;
                 dblTemp = (double*)malloc(iSize * sizeof(double));
                 C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne);
 
-                cNorm = 'N';
-                C2F(dgetrs)(&cNorm, &_iCols1, &_iCols2, pAf, &_iCols1, pIpiv, dblTemp, &_iCols1, &iInfo);
-                cNorm = 'F';
-                C2F(dlacpy)(&cNorm, &_iCols1, &_iCols2, dblTemp, &_iCols1, _pdblRealOut, &_iCols1);
+                //solve triangular system
+                C2F(dtrsm)("L", &cUpLo,"N", "N", &_iRows2, &_iCols2, &dblOne, _pdblReal1, &_iRows1, dblTemp, &_iRows2);
+                C2F(dcopy)(&iSize, dblTemp, &iOne, _pdblRealOut, &iOne);
                 iExit = 1;
 
                 free(dblTemp);
             }
         }
+        else
+        {
+            iSize = _iRows1 * _iCols1;
+            C2F(dcopy)(&iSize, _pdblReal1, &iOne, pAf, &iOne);
+            C2F(dgetrf)(&_iCols1, &_iCols1, pAf, &_iCols1, pIpiv, &iInfo);
+            if (iInfo == 0)
+            {
+                dblAnorm = C2F(dlange)("1", &_iRows1, &_iCols1, _pdblReal1, &_iRows1, pDwork);
+                C2F(dgecon)("1", &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
+                if (dblRcond > RCONDthresh)
+                {
+                    // _pdblReal2 will be overwritten by dgetrs
+                    iSize = _iRows2 * _iCols2;
+                    dblTemp = (double*)malloc(iSize * sizeof(double));
+                    C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne);
+
+                    C2F(dgetrs)("N", &_iCols1, &_iCols2, pAf, &_iCols1, pIpiv, dblTemp, &_iCols1, &iInfo);
+                    C2F(dcopy)(&iSize, dblTemp, &iOne, _pdblRealOut, &iOne);
+                    iExit = 1;
+
+                    free(dblTemp);
+                }
+            }
+        }
 
         if (iExit == 0)
         {
@@ -646,7 +671,7 @@ int iLeftDivisionOfRealMatrix(
         iMax = Max(_iRows1, _iCols1);
         C2F(dlacpy)(&cNorm, &_iRows1, &_iCols2, _pdblReal2, &_iRows1, pXb, &iMax);
         memset(pJpvt, 0x00, sizeof(int) * _iCols1);
-        // _pdblReal1 will be overwrite by dgelsy1
+        // _pdblReal1 will be overwritten by dgelsy1
         iSize = _iRows1 * _iCols1;
         dblTemp = (double*)malloc(iSize * sizeof(double));
         C2F(dcopy)(&iSize, _pdblReal1, &iOne, dblTemp, &iOne);
@@ -689,6 +714,8 @@ int iLeftDivisionOfComplexMatrix(
     int iIndex  = 0;
     char cNorm  = 0;
     int iExit   = 0;
+    int iLower = 0;
+    int iUpper = 0;
 
     /*temporary variables*/
     int iWorkMin    = 0;
@@ -699,6 +726,7 @@ int iLeftDivisionOfComplexMatrix(
     double dblEps       = 0;
     double RCONDthresh  = 0;
     double dblAnorm     = 0;
+    doublecomplex dblCplxOne = {1,0};
 
     doublecomplex *pAf      = NULL;
     doublecomplex *pXb      = NULL;
@@ -727,18 +755,20 @@ int iLeftDivisionOfComplexMatrix(
     pDwork      = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
     dblEps      = nc_eps();
     RCONDthresh = 10 * dblEps;
-    dblAnorm    = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
 
     if (_iRows1 == _iCols1)
     {
-        C2F(zgetrf)(&_iCols1, &_iCols1, poVar1, &_iCols1, pIpiv, &iInfo);
-        if (iInfo == 0)
+        matrixIsTriangular(_pdblReal1, _pdblImg1, _iRows1, _iCols1, &iUpper, &iLower);
+        if (iUpper || iLower)
         {
-            C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
+            //if matrix is triangular
+            char cUpLoType[4] = {'N','U','L','U'};
+            char cUpLo = cUpLoType[iUpper + 2*iLower];
+            C2F(ztrcon)("1", &cUpLo, "N", &_iRows1, poVar1, &_iRows1, &dblRcond, pDwork, pRwork, &iInfo);
             if (dblRcond > RCONDthresh)
             {
-                cNorm    = 'N';
-                C2F(zgetrs)(&cNorm, &_iCols1, &_iCols2, poVar1, &_iCols1, pIpiv, poVar2, &_iCols1, &iInfo);
+                //solve triangular system
+                C2F(ztrsm)("L", &cUpLo,"N", "N", &_iRows2, &_iCols2, &dblCplxOne, poVar1, &_iRows1, poVar2, &_iRows2);
                 vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
                 iExit = 1;
             }
@@ -749,6 +779,27 @@ int iLeftDivisionOfComplexMatrix(
                 *_pdblRcond = dblRcond;
             }
         }
+        else
+        {
+            C2F(zgetrf)(&_iCols1, &_iCols1, poVar1, &_iCols1, pIpiv, &iInfo);
+            if (iInfo == 0)
+            {
+                dblAnorm = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
+                C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
+                if (dblRcond > RCONDthresh)
+                {
+                    C2F(zgetrs)("N", &_iCols1, &_iCols2, poVar1, &_iCols1, pIpiv, poVar2, &_iCols1, &iInfo);
+                    vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
+                    iExit = 1;
+                }
+                else
+                {
+                    //how to extract that ? Oo
+                    iReturn = -1;
+                    *_pdblRcond = dblRcond;
+                }
+            }
+        }
     }
 
     if (iExit == 0)
@@ -798,3 +849,44 @@ int iLeftDivisionOfComplexMatrix(
     free(pDwork);
     return 0;
 }
+
+void matrixIsTriangular(double *_pdblReal, double *_pdblImg, int _iRows, int _iCols, int *bUpper, int *bLower)
+{
+    double *pdblArray[2] = {_pdblReal, _pdblImg};
+    double *pdbl;
+
+    *bUpper = 1;
+    *bLower = 1;
+
+    for (int i=0; i<2 && pdblArray[i] != NULL; i++)
+    {
+        int j;
+        int iDim;
+        int iOne = 1;
+        //upper triangular ?
+        if (*bUpper)
+        {
+            pdbl = pdblArray[i] + 1;
+            iDim = _iRows-1;
+            for (j = 0; j < _iCols; j++, iDim--)
+            {
+                // compute L1 norm with BLAS appears to be faster than testing
+                // successive nullity of individual terms
+                if (C2F(dasum)(&iDim, pdbl, &iOne) > 0.0) break;
+                pdbl += _iRows+1;
+            }
+            *bUpper = (j==_iCols);
+        }
+        //lower triangular ?
+        if (*bLower)
+        {
+            pdbl = pdblArray[i] + _iRows;
+            for (j = 1; j < _iCols; j++)
+            {
+                if (C2F(dasum)(&j, pdbl, &iOne) > 0.0) break;
+                pdbl += _iRows;
+            }
+            *bLower = (j==_iCols);
+        }
+    }
+}
diff --git a/scilab/modules/ast/tests/nonreg_tests/bug_15701.tst b/scilab/modules/ast/tests/nonreg_tests/bug_15701.tst
new file mode 100644 (file)
index 0000000..b67dfe1
--- /dev/null
@@ -0,0 +1,57 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2019 - Stéphane MOTTELET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+//
+// <-- Non-regression test for bug 15701 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/15701
+//
+// <-- Short Description -->
+// A\B is not faster when A is square and triangular
+
+N = 3000;
+// matrix A with good condition number in order to prevent
+// least squares solution
+A=tril(rand(N,N))+10*eye(N,N);
+B=rand(N,1);
+
+// general case
+A(1,N)=%eps;
+tic;
+x1=A\B;
+t1=toc();
+
+// triangular case
+A(1,N)=0;
+tic;
+x2=A\B;
+t2=toc();
+
+assert_checkalmostequal(x1,x2);
+assert_checktrue(t1/t2 > 10);
+
+// complex case
+A=A+%i*tril(rand(N,N));
+B=B+%i*rand(N,1);
+
+// general case
+A(1,N)=%eps;
+tic;
+x1=A\B;
+t1=toc();
+
+// triangular case
+A(1,N)=0;
+tic;
+x2=A\B;
+t2=toc();
+
+assert_checkalmostequal(x1,x2);
+assert_checktrue(t1/t2 > 10);
\ No newline at end of file
index 74890fd..eb2aa16 100644 (file)
@@ -73,6 +73,8 @@ extern int C2F(zscal) (int *_iSize, doublecomplex * _pdblVal, doublecomplex * _p
 extern int C2F(dset) (int *_iSize, double *_pdblVal, double *_pdblDest, int *_iInc);
 extern double C2F(dlange) (char *_pstNorm, int *_piM, int *_piN, double *_pdblA, int *_piLdA, double *_pdblWork);
 extern int C2F(dlacpy) (char *_pstUplo, int *piM, int *_piN, double *_pdblA, int *_piLdA, double *_pdblB, int *_piLdB);
+extern int C2F(dtrcon) (char *_pstNORM, char*uplo, char *diag, int *_piN, double *_pdblA, int *_piLDA, double *_pdblRCOND, double *_pdblWORK,
+                        int *_piIWORK, int *_piINFO);
 extern int C2F(dgecon) (char *_pstNORM, int *_piN, double *_pdblA, int *_piLDA, double *_pdblANORM, double *_pdblRCOND, double *_pdblWORK,
                         int *_piIWORK, int *_piINFO);
 extern int C2F(dgelsy1) (int *_piM, int *_piN, int *_piNRHS, double *_pdblA, int *_piLDA, double *_pdblB, int *_piLDB, int *_piJPVT,
@@ -81,6 +83,8 @@ extern double C2F(zlange) (char *_pstNORM, int *_piM, int *_piN, double *_pdblA,
 extern int C2F(zlacpy) (char *_pstUPLO, int *_piM, int *_piN, double *_pdblA, int *_piLDA, double *_pdblB, int *_piLDB);
 extern void C2F(zgecon) (char *_pstNORM, int *_piN, doublecomplex * _pdblA, int *_piLDA, double *_pdblANORM, double *_pdblRNORM,
                          doublecomplex * _pdblWORK, double *_pdblRWORD, int *_piINFO);
+extern int C2F(ztrcon) (char *_pstNORM, char*uplo, char *diag, int *_piN, doublecomplex *_pdblA, int *_piLDA, double *_pdblRCOND, double *_pdblWORK,
+                                                 double *_pdblRWORD, int *_piINFO);
 extern int C2F(zgelsy1) (int *_piM, int *_piN, int *_piNRHS, doublecomplex * pdblA, int *_piLDA, doublecomplex * _pdblB, int *_piLDB, int *_piJPVT,
                          double *_pdblRCOND, int *_piRANK, doublecomplex * _pdblWORK, int *_piLWORK, double *_pdblRWORK, int *_piINFO);
 extern double C2F(ddot) (int *_ipSize, double *_pdblVal1, int *_piInc1, double *_pdblVal2, int *_piInc2);