complex left division corrected. 62/11562/5
Cedric Delamarre [Wed, 22 May 2013 15:28:01 +0000 (17:28 +0200)]
a=rand(4,4);b=rand(5,4);ac=a+%i*rand(4,4);
ac\b'
a\%i

test_run("core","LU",["no_check_error_output" ]);
test_run("core","QR",["no_check_error_output" ]);

Change-Id: I59dba8ab0186c05577f4178e088adcd645661f07

scilab/modules/operations/src/c/matrix_division.c
scilab/modules/operations/src/cpp/types_ldivide.cpp

index 4279731..815b4c3 100644 (file)
@@ -698,50 +698,36 @@ int       iLeftDivisionOfComplexMatrix(
 
     double *pRwork                     = NULL;
 
-    int *pRank                         = NULL;
+    int iRank                          = 0;
     int *pIpiv                         = NULL;
     int *pJpvt                         = NULL;
 
     iWorkMin   = Max(2 * _iCols1, Min(_iRows1, _iCols1) + Max(2 * Min(_iRows1, _iCols1), Max(_iCols1, Min(_iRows1, _iCols1) + _iCols2)));
 
     /* Array allocations*/
-    poVar1             = oGetDoubleComplexFromPointer(_pdblReal1,      _pdblImg1,              _iRows1 * _iCols1);
-    poVar2             = oGetDoubleComplexFromPointer(_pdblReal2,      _pdblImg2,              _iRows2 * _iCols2);
-    poOut              = oGetDoubleComplexFromPointer(_pdblRealOut, _pdblImgOut,       _iRowsOut * _iColsOut);
+    poVar1         = oGetDoubleComplexFromPointer(_pdblReal1, _pdblImg1, _iRows1 * _iCols1);
+    poVar2      = oGetDoubleComplexFromPointer(_pdblReal2, _pdblImg2, _iRows2 * _iCols2);
 
-    pAf                        = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows1 * _iCols1);
-    pXb                        = (doublecomplex*)malloc(sizeof(doublecomplex) * Max(_iRows1, _iCols1) * _iColsOut);
+    pIpiv       = (int*)malloc(sizeof(int) * _iCols1);
+    pJpvt       = (int*)malloc(sizeof(int) * _iCols1);
+    pRwork      = (double*)malloc(sizeof(double) * _iCols1 * 2);
 
-    pRank              = (int*)malloc(sizeof(int));
-    pIpiv              = (int*)malloc(sizeof(int) * _iCols1);
-    pJpvt              = (int*)malloc(sizeof(int) * _iCols1);
-    pRwork             = (double*)malloc(sizeof(double) * _iCols1 * 2);
-
-
-    //C'est du grand nawak ca, on reserve toute la stack ! Oo
-
-    cNorm              = '1';
-    pDwork             = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
-    dblEps             = F2C(dlamch)("eps", 1L);
-    dblAnorm   = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
+    cNorm       = '1';
+    pDwork      = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
+    dblEps      = F2C(dlamch)("eps", 1L);
+    dblAnorm    = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
 
     if (_iRows1 == _iCols1)
     {
-        cNorm          = 'F';
-        C2F(zlacpy)(&cNorm, &_iCols1, &_iCols1,        (double*)poVar1, &_iCols1, (double*)pAf, &_iCols1);
-        C2F(zlacpy)(&cNorm, &_iCols1, &_iCols2,        (double*)poVar2, &_iCols1, (double*)pXb, &_iCols1);
-        C2F(zgetrf)(&_iCols1, &_iCols1, (double*)pAf, &_iCols1, pIpiv, &iInfo);
+        C2F(zgetrf)(&_iCols1, &_iCols1, (double*)poVar1, &_iCols1, pIpiv, &iInfo);
         if (iInfo == 0)
         {
-            cNorm = '1';
-            C2F(zgecon)(&cNorm, &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
+            C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
             if (dblRcond > dsqrts(dblEps))
             {
                 cNorm  = 'N';
-                C2F(zgetrs)(&cNorm, &_iCols1, &_iCols2, (double*)pAf, &_iCols1, pIpiv, (double*)pXb, &_iCols1, &iInfo);
-                cNorm  = 'F';
-                C2F(zlacpy)(&cNorm, &_iCols1, &_iCols2, (double*)pXb, &_iCols1, (double*)poOut, &_iCols1);
-                vGetPointerFromDoubleComplex(poOut, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
+                C2F(zgetrs)(&cNorm, &_iCols1, &_iCols2, (double*)poVar1, &_iCols1, pIpiv, (double*)poVar2, &_iCols1, &iInfo);
+                vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
                 iExit = 1;
             }
             else
@@ -756,36 +742,32 @@ int       iLeftDivisionOfComplexMatrix(
     if (iExit == 0)
     {
         dblRcond = dsqrts(dblEps);
-        cNorm = 'F';
         iMax = Max(_iRows1, _iCols1);
-        C2F(zlacpy)(&cNorm, &_iRows1, &_iCols2, (double*)poVar2, &_iRows1, (double*)pXb, &iMax);
         memset(pJpvt, 0x00, sizeof(int) * _iCols1);
-        C2F(zgelsy1)(  &_iRows1, &_iCols1, &_iCols2, poVar1, &_iRows1, pXb, &iMax,
-                        pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, pRwork, &iInfo);
+        pXb = (doublecomplex*)malloc(sizeof(doublecomplex) * iMax * _iColsOut);
+        cNorm = 'F';
+        C2F(zlacpy)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)pXb, &iMax);
+        // pXb : in input pXb is of size rows1 x col2
+        //       in output pXp is of size col1 x col2
+        C2F(zgelsy1)(&_iRows1, &_iCols1, &_iCols2, poVar1, &_iRows1, pXb, &iMax,
+                     pJpvt, &dblRcond, &iRank, pDwork, &iWorkMin, pRwork, &iInfo);
 
         if (iInfo == 0)
         {
-            if ( _iRows1 != _iCols1 && pRank[0] < Min(_iRows1, _iCols1))
+            if ( _iRows1 != _iCols1 && iRank < Min(_iRows1, _iCols1))
             {
                 //how to extract that ? Oo
                 iReturn = -2;
-                *_pdblRcond = pRank[0];
+                *_pdblRcond = (double)iRank;
             }
 
-            cNorm = 'F';
-            C2F(zlacpy)(&cNorm, &_iCols1, &_iCols2, (double*)pXb, &iMax, (double*)poOut, &_iCols1);
-            vGetPointerFromDoubleComplex(poOut, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
+            vGetPointerFromDoubleComplex(pXb, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
         }
+        free(pXb);
     }
 
-
     vFreeDoubleComplexFromPointer(poVar1);
     vFreeDoubleComplexFromPointer(poVar2);
-    vFreeDoubleComplexFromPointer(poOut);
-
-    free(pAf);
-    free(pXb);
-    free(pRank);
     free(pIpiv);
     free(pJpvt);
     free(pRwork);
index 5ace5a1..8588fe7 100644 (file)
@@ -173,7 +173,7 @@ int LDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubl
         iErr = iLeftDivisionOfComplexMatrix(
                    _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
                    _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
-                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
+                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), (*_pDoubleOut)->getRows(), (*_pDoubleOut)->getCols(), &dblRcond);
     }
     else
     {
@@ -181,7 +181,7 @@ int LDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubl
         iErr = iLeftDivisionOfRealMatrix(
                    _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
                    _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
-                   (*_pDoubleOut)->getReal(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
+                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getRows(), (*_pDoubleOut)->getCols(), &dblRcond);
     }
 
     return iErr;