* Bug 15701 fixed: now A\B is faster when A is square and triangular
[scilab.git] / scilab / modules / ast / src / c / operations / matrix_division.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2008-2008 - INRIA - Antoine ELIAS <antoine.elias@scilab.org>
4  * Copyright (C) 2019 - St├ęphane MOTTELET
5  *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  *
15  */
16
17 #include "configvariable_interface.h"
18 #include "matrix_division.h"
19 #include <string.h>
20
21 /*iRightDivisionComplexMatrixByComplexMatrix*/
22 int iRightDivisionComplexMatrixByComplexMatrix(
23     double *_pdblReal1,     double *_pdblImg1,      int _iInc1,
24     double *_pdblReal2,     double *_pdblImg2,      int _iInc2,
25     double *_pdblRealOut,   double *_pdblImgOut,    int _iIncOut,   int _iSize)
26 {
27     int iErr        = 0;
28     int iIndex      = 0; //Main loop index
29     int iIndex1     = 0; //Loop index on left operand
30     int iIndex2     = 0; //Loop index on right operand
31     int iIndexOut   = 0; //Lopp index on result matrix
32
33     if (_iInc2 == 0)
34     {
35         if ((getieee() == 0) && (dabss(_pdblReal2[iIndex2]) + dabss(_pdblImg2[iIndex2]) == 0))
36         {
37             return 3;
38         }
39     }
40
41     for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
42     {
43         iErr = iRightDivisionComplexByComplex(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
44         iIndexOut    += _iIncOut;
45         iIndex1      += _iInc1;
46         iIndex2      += _iInc2;
47     }
48
49     return iErr;
50 }
51
52 /*iRightDivisionComplexByComplex*/
53 int iRightDivisionComplexByComplex(
54     double _dblReal1, double _dblImg1,
55     double _dblReal2, double _dblImg2,
56     double *_pdblRealOut, double *_pdblImgOut)
57 {
58     int iErr = 0;
59     if (_dblImg2 == 0)
60     {
61         if (_dblReal2 == 0)
62         {
63             //got NaN + i NaN
64             iErr = 10;
65             *_pdblRealOut   = _dblImg2 / _dblReal2;
66             *_pdblImgOut    = *_pdblRealOut;
67         }
68         else
69         {
70             *_pdblRealOut   = _dblReal1 / _dblReal2;
71             *_pdblImgOut    = _dblImg1 / _dblReal2;
72         }
73     }
74     else if (_dblReal2 == 0)
75     {
76         *_pdblRealOut   = _dblImg1 / _dblImg2;
77         *_pdblImgOut    = (-_dblReal1) / _dblImg2;
78     }
79     else
80     {
81         //Generic division algorithm
82
83         if (dabss(_dblReal2) >= dabss(_dblImg2))
84         {
85             double dblRatio2    = _dblImg2 / _dblReal2;
86             double dblSum       = _dblReal2 + dblRatio2 * _dblImg2;
87             *_pdblRealOut       = (_dblReal1 + _dblImg1 * dblRatio2) / dblSum;
88             *_pdblImgOut        = (_dblImg1 - _dblReal1 * dblRatio2) / dblSum;
89         }
90         else
91         {
92             double dblRatio2    = _dblReal2 / _dblImg2;
93             double dblSum       = _dblImg2 +  dblRatio2 * _dblReal2;
94             *_pdblRealOut       = (_dblReal1 * dblRatio2 + _dblImg1) / dblSum;
95             *_pdblImgOut        = (_dblImg1 * dblRatio2 - _dblReal1) / dblSum;
96         }
97     }
98     return iErr;
99 }
100
101 /*iRightDivisionRealMatrixByComplexMatrix*/
102 int iRightDivisionRealMatrixByComplexMatrix(
103     double *_pdblReal1,     int _iInc1,
104     double *_pdblReal2,     double *_pdblImg2,      int _iInc2,
105     double *_pdblRealOut,   double *_pdblImgOut,    int _iIncOut, int _iSize)
106 {
107     int iErr        = 0;
108     int iIndex      = 0; //Main loop index
109     int iIndex1     = 0; //Loop index on left operand
110     int iIndex2     = 0; //Loop index on right operand
111     int iIndexOut   = 0; //Lopp index on result matrix
112
113     if (_iInc2 == 0)
114     {
115         if ((getieee() == 0) && (dabss(_pdblReal2[iIndex2]) + dabss(_pdblImg2[iIndex2]) == 0))
116         {
117             return 3;
118         }
119     }
120
121     for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
122     {
123         iErr = iRightDivisionRealByComplex(_pdblReal1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
124         iIndexOut    += _iIncOut;
125         iIndex1      += _iInc1;
126         iIndex2      += _iInc2;
127     }
128
129     return iErr;
130 }
131
132 /*iRightDivisionRealByComplex*/
133 int iRightDivisionRealByComplex(
134     double _dblReal1,
135     double _dblReal2, double _dblImg2,
136     double *_pdblRealOut, double *_pdblImgOut)
137 {
138     int iErr = 0;
139     if (_dblImg2 == 0)
140     {
141         *_pdblRealOut   = _dblReal1 / _dblReal2;
142         *_pdblImgOut    = 0;
143     }
144     else if (_dblReal2 == 0)
145     {
146         *_pdblRealOut   = 0;
147         *_pdblImgOut    = -_dblReal1 / _dblImg2;
148     }
149     else
150     {
151         double dblAbsSum = dabss(_dblReal2) + dabss(_dblImg2);
152
153         if (dblAbsSum == 0)
154         {
155             iErr = 10;
156             *_pdblRealOut   = _dblReal1 / dblAbsSum;
157             *_pdblImgOut    = 0;
158         }
159         else
160         {
161             double dblReal1Sum  = _dblReal1 / dblAbsSum;
162             double dblReal2Sum  = _dblReal2 / dblAbsSum;
163             double dblImg2Sum   = _dblImg2 / dblAbsSum;
164             double dblSum       = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2);
165             *_pdblRealOut       = (dblReal1Sum * dblReal2Sum) / dblSum;
166             *_pdblImgOut        = (-dblReal1Sum * dblImg2Sum) / dblSum;
167         }
168     }
169     return iErr;
170 }
171
172 /*iRightDivisionComplexMatrixByRealMatrix*/
173 int iRightDivisionComplexMatrixByRealMatrix(
174     double *_pdblReal1,        double *_pdblImg1,        int _iInc1,
175     double *_pdblReal2,                                int _iInc2,
176     double *_pdblRealOut,    double *_pdblImgOut,    int _iIncOut, int _iSize)
177 {
178     int iErr        = 0;
179     int iIndex      = 0; //Main loop index
180     int iIndex1     = 0; //Loop index on left operand
181     int iIndex2     = 0; //Loop index on right operand
182     int iIndexOut   = 0; //Lopp index on result matrix
183     for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
184     {
185         iErr = iRightDivisionComplexByReal(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
186         iIndexOut   += _iIncOut;
187         iIndex1     += _iInc1;
188         iIndex2     += _iInc2;
189     }
190     return iErr;
191 }
192
193 /*iRightDivisionComplexByReal*/
194 int iRightDivisionComplexByReal(
195     double _dblReal1, double _dblImg1,
196     double _dblReal2,
197     double *_pdblRealOut, double *_pdblImgOut)
198 {
199     int iErr = 0;
200     if (_dblReal2 == 0)
201     {
202         if (getieee() == 0)
203         {
204             iErr = 3;
205             return iErr;
206         }
207         else if (getieee() == 1)
208         {
209             //Warning
210             iErr = 4;
211         }
212     }
213
214     *_pdblRealOut    = _dblReal1 / _dblReal2;
215     *_pdblImgOut    = _dblImg1 / _dblReal2;
216
217     return iErr;
218 }
219
220 /*iRightDivisionRealMatrixByRealMatrix*/
221 int iRightDivisionRealMatrixByRealMatrix(
222     double *_pdblReal1, int _iInc1,
223     double *_pdblReal2, int _iInc2,
224     double *_pdblRealOut, int _iIncOut, int _iSize)
225 {
226     int iIndex      = 0; //Main loop index
227     int iIndex1     = 0; //Loop index on left operand
228     int iIndex2     = 0; //Loop index on right operand
229     int iIndexOut   = 0; //Lopp index on result matrix
230     int iErr        = 0;
231
232     for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
233     {
234         if (_pdblReal2[iIndex2] == 0)
235         {
236             if (getieee() == 0)
237             {
238                 iErr = 3;
239                 return iErr;
240             }
241             else if (getieee() == 1)
242             {
243                 iErr = 4;
244             }
245         }
246
247         _pdblRealOut[iIndexOut] = _pdblReal1[iIndex1] / _pdblReal2[iIndex2];
248         iIndexOut              += _iIncOut;
249         iIndex1                += _iInc1;
250         iIndex2                += _iInc2;
251     }
252     return iErr;
253 }
254
255 int iRightDivisionOfRealMatrix(
256     double *_pdblReal1,     int _iRows1,    int _iCols1,
257     double *_pdblReal2,     int _iRows2,    int _iCols2,
258     double *_pdblRealOut,   int _iRowsOut,  int _iColsOut,  double* _pdblRcond)
259 {
260     int iReturn = 0;
261
262     int iIndex  = 0;
263     char cNorm  = 0;
264     int iExit   = 0;
265
266     /*temporary variables*/
267     int iWorkMin    = 0;
268     int iInfo       = 0;
269     int iMax        = 0;
270     double dblRcond = 0;
271
272     double dblEps       = 0;
273     double RCONDthresh  = 0;
274     double dblAnorm     = 0;
275
276     double *pAf     = NULL;
277     double *pAt     = NULL;
278     double *pBt     = NULL;
279     double *pDwork  = NULL;
280
281     int *pRank    = NULL;
282     int *pIpiv    = NULL;
283     int *pJpvt    = NULL;
284     int *pIwork   = NULL;
285
286     iWorkMin    = Max(4 * _iCols2, Max(Min(_iRows2, _iCols2) + 3 * _iRows2 + 1, 2 * Min(_iRows2, _iCols2) + _iRows1));
287
288     /* Array allocations*/
289     pAf         = (double*)malloc(sizeof(double) * _iCols2 * _iRows2);
290     pAt         = (double*)malloc(sizeof(double) * _iCols2 * _iRows2);
291     pBt         = (double*)malloc(sizeof(double) * Max(_iRows2, _iCols2) * _iRows1);
292
293     pRank       = (int*)malloc(sizeof(int));
294     pIpiv       = (int*)malloc(sizeof(int) * _iCols2);
295     pJpvt       = (int*)malloc(sizeof(int) * _iRows2);
296     pIwork      = (int*)malloc(sizeof(int) * _iCols2);
297
298
299     //C'est du grand nawak ca, on reserve toute la stack ! Oo
300
301     cNorm       = '1';
302     pDwork      = (double*)malloc(sizeof(double) * iWorkMin);
303     dblEps      = nc_eps();
304     RCONDthresh = 10 * dblEps;
305     dblAnorm    = C2F(dlange)(&cNorm, &_iRows2, &_iCols2, _pdblReal2, &_iRows2, pDwork);
306
307     //tranpose A and B
308
309     vTransposeRealMatrix(_pdblReal2, _iRows2, _iCols2, pAt);
310
311     {
312         int i, j, ij, ji;
313         for (j = 0 ; j < _iRows1 ; j++)
314         {
315             for (i = 0 ; i < _iCols2 ; i++)
316             {
317                 ij = i + j * Max(_iRows2, _iCols2);
318                 ji = j + i * _iRows1;
319                 pBt[ij]    = _pdblReal1[ji];
320             }//for(j = 0 ; j < _iRows1 ; j++)
321         }//for(i = 0 ; i < _iCols2 ; i++)
322     }//bloc esthetique
323
324     if (_iRows2 == _iCols2)
325     {
326         cNorm = 'F';
327         C2F(dlacpy)(&cNorm, &_iCols2, &_iCols2, pAt, &_iCols2, pAf, &_iCols2);
328         C2F(dgetrf)(&_iCols2, &_iCols2, pAf, &_iCols2, pIpiv, &iInfo);
329         if (iInfo == 0)
330         {
331             cNorm = '1';
332             C2F(dgecon)(&cNorm, &_iCols2, pAf, &_iCols2, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
333             if (dblRcond > RCONDthresh)
334             {
335                 cNorm = 'N';
336                 C2F(dgetrs)(&cNorm, &_iCols2, &_iRows1, pAf, &_iCols2, pIpiv, pBt, &_iCols2, &iInfo);
337                 vTransposeRealMatrix(pBt, _iCols2, _iRows1, _pdblRealOut);
338                 iExit = 1;
339             }
340         }
341
342         if (iExit == 0)
343         {
344             //how to extract that ? Oo
345             *_pdblRcond = dblRcond;
346             iReturn = -1;
347         }
348     }
349
350     if (iExit == 0)
351     {
352         dblRcond = RCONDthresh;
353         cNorm = 'F';
354         iMax = Max(_iRows2, _iCols2);
355         memset(pJpvt, 0x00, sizeof(int) * _iRows2);
356         iInfo = 1;
357         C2F(dgelsy1)(&_iCols2, &_iRows2, &_iRows1, pAt, &_iCols2, pBt, &iMax,
358                      pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo);
359
360         if (iInfo == 0)
361         {
362             if ( _iRows2 != _iCols2 && pRank[0] < Min(_iRows2, _iCols2))
363             {
364                 //how to extract that ? Oo
365                 iReturn = -2;
366                 *_pdblRcond = pRank[0];
367             }
368
369             //    TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2);
370
371             //Mega caca de la mort qui tue des ours a mains nues
372             //mais je ne sais pas comment le rendre "beau" :(
373             {
374                 int i, j, ij, ji;
375                 for (j = 0 ; j < _iRows2 ; j++)
376                 {
377                     for (i = 0 ; i < _iRows1 ; i++)
378                     {
379                         ij = i + j * _iRows1;
380                         ji = j + i * Max(_iRows2, _iCols2);
381                         _pdblRealOut[ij]    = pBt[ji];
382                     }//for(i = 0 ; i < _iRows2 ; i++)
383                 }//for(j = 0 ; j < _iRows1 ; j++)
384             }//bloc esthetique
385         }//if(iInfo == 0)
386     }//if(bExit == 0)
387
388     free(pAf);
389     free(pAt);
390     free(pBt);
391     free(pRank);
392     free(pIpiv);
393     free(pJpvt);
394     free(pIwork);
395     free(pDwork);
396     return iReturn;
397 }
398
399
400 int iRightDivisionOfComplexMatrix(
401     double *_pdblReal1,     double *_pdblImg1,      int _iRows1,    int _iCols1,
402     double *_pdblReal2,     double *_pdblImg2,      int _iRows2,    int _iCols2,
403     double *_pdblRealOut,   double *_pdblImgOut,    int _iRowsOut,  int _iColsOut,  double *_pdblRcond)
404 {
405     int iReturn     = 0;
406     int iIndex1     = 0;
407     int iIndex2     = 0;
408     char cNorm      = 0;
409     int iExit       = 0;
410
411     /*temporary variables*/
412     int iWorkMin    = 0;
413     int iInfo       = 0;
414     int iMax        = 0;
415     double dblRcond = 0;
416
417     double dblEps       = 0;
418     double RCONDthresh  = 0;
419     double dblAnorm     = 0;
420
421     doublecomplex *poVar1   = NULL;
422     doublecomplex *poVar2   = NULL;
423     doublecomplex *poOut    = NULL;
424     doublecomplex *poAf     = NULL;
425     doublecomplex *poAt     = NULL;
426     doublecomplex *poBt     = NULL;
427     doublecomplex *poDwork  = NULL;
428
429     int *pRank    = NULL;
430     int *pIpiv    = NULL;
431     int *pJpvt    = NULL;
432     double *pRwork    = NULL;
433
434     iWorkMin    = Max(2 * _iCols2, Min(_iRows2, _iCols2) + Max(2 * Min(_iRows2, _iCols2), Max(_iRows2 + 1, Min(_iRows2, _iCols2) + _iRows1)));
435
436     /* Array allocations*/
437     poVar1      = oGetDoubleComplexFromPointer(_pdblReal1,        _pdblImg1,        _iRows1 * _iCols1);
438     poVar2      = oGetDoubleComplexFromPointer(_pdblReal2,        _pdblImg2,        _iRows2 * _iCols2);
439     poOut       = oGetDoubleComplexFromPointer(_pdblRealOut,    _pdblImgOut,    _iRowsOut * _iColsOut);
440
441     poAf        = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2);
442     poAt        = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2);
443     poBt        = (doublecomplex*)malloc(sizeof(doublecomplex) * Max(_iRows2, _iCols2) * _iRows1);
444     poDwork     = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
445
446     pRank       = (int*)malloc(sizeof(int));
447     pIpiv       = (int*)malloc(sizeof(int) * _iCols2);
448     pJpvt       = (int*)malloc(sizeof(int) * _iRows2);
449     pRwork      = (double*)malloc(sizeof(double) * 2 * _iRows2);
450
451     dblEps      = nc_eps();
452     RCONDthresh = 10 * dblEps;
453     cNorm       = '1';
454     dblAnorm    = C2F(zlange)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)poDwork);
455
456     //tranpose A and B
457
458     vTransposeDoubleComplexMatrix(poVar2, _iRows2, _iCols2, poAt, 1);
459
460     {
461         int i, j, ij, ji;
462         for (j = 0 ; j < _iRows1 ; j++)
463         {
464             for (i = 0 ; i < _iCols2 ; i++)
465             {
466                 ij = i + j * Max(_iRows2, _iCols2);
467                 ji = j + i * _iRows1;
468                 poBt[ij].r    = poVar1[ji].r;
469                 //Conjugate
470                 poBt[ij].i    = -poVar1[ji].i;
471             }//for(j = 0 ; j < _iRows1 ; j++)
472         }//for(i = 0 ; i < _iCols2 ; i++)
473     }//bloc esthetique
474
475
476     if (_iRows2 == _iCols2)
477     {
478         cNorm = 'F';
479         C2F(zlacpy)(&cNorm, &_iCols2, &_iCols2, (double*)poAt, &_iCols2, (double*)poAf, &_iCols2);
480         C2F(zgetrf)(&_iCols2, &_iCols2, poAf, &_iCols2, pIpiv, &iInfo);
481         if (iInfo == 0)
482         {
483             cNorm = '1';
484             C2F(zgecon)(&cNorm, &_iCols2, poAf, &_iCols2, &dblAnorm, &dblRcond, poDwork, pRwork, &iInfo);
485             if (dblRcond > RCONDthresh)
486             {
487                 cNorm = 'N';
488                 C2F(zgetrs)(&cNorm, &_iCols2, &_iRows1, poAf, &_iCols2, pIpiv, poBt, &_iCols2, &iInfo);
489                 vTransposeDoubleComplexMatrix(poBt, _iCols2, _iRows1, poOut, 1);
490                 vGetPointerFromDoubleComplex(poOut, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
491                 iExit = 1;
492             }
493         }
494
495         if (iExit == 0)
496         {
497             //how to extract that ? Oo
498             *_pdblRcond = dblRcond;
499             iReturn = -1;
500         }
501     }
502
503     if (iExit == 0)
504     {
505         dblRcond = RCONDthresh;
506         cNorm = 'F';
507         iMax = Max(_iRows2, _iCols2);
508         memset(pJpvt, 0x00, sizeof(int) * _iRows2);
509         iInfo = 1;
510         C2F(zgelsy1)(&_iCols2, &_iRows2, &_iRows1, poAt, &_iCols2, poBt, &iMax,
511                      pJpvt, &dblRcond, pRank, poDwork, &iWorkMin, pRwork, &iInfo);
512
513         if (iInfo == 0)
514         {
515             if ( _iRows2 != _iCols2 && pRank[0] < Min(_iRows2, _iCols2))
516             {
517                 //how to extract that ? Oo
518                 iReturn = -2;
519                 *_pdblRcond = pRank[0];
520             }
521
522             //    TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2);
523
524             //Mega caca de la mort qui tue des ours a mains nues
525             //mais je ne sais pas comment le rendre "beau" :(
526             {
527                 int i, j, ij, ji;
528                 for (j = 0 ; j < _iRows2 ; j++)
529                 {
530                     for (i = 0 ; i < _iRows1 ; i++)
531                     {
532                         ij = i + j * _iRows1;
533                         ji = j + i * Max(_iRows2, _iCols2);
534                         _pdblRealOut[ij]    = poBt[ji].r;
535                         //Conjugate
536                         _pdblImgOut[ij]        = -poBt[ji].i;
537                     }//for(i = 0 ; i < _iRows2 ; i++)
538                 }//for(j = 0 ; j < _iRows1 ; j++)
539             }//bloc esthetique
540         }//if(iInfo == 0)
541     }//if(iExit == 0)
542
543
544     vFreeDoubleComplexFromPointer(poVar1);
545     vFreeDoubleComplexFromPointer(poVar2);
546     vFreeDoubleComplexFromPointer(poOut);
547
548     free(poAf);
549     free(poAt);
550     free(poBt);
551     free(pRank);
552     free(pIpiv);
553     free(pJpvt);
554     free(pRwork);
555     free(poDwork);
556     return 0;
557 }
558
559 /*Matrix left division*/
560 int iLeftDivisionOfRealMatrix(
561     double *_pdblReal1,     int _iRows1,    int _iCols1,
562     double *_pdblReal2,     int _iRows2,    int _iCols2,
563     double *_pdblRealOut,   int _iRowsOut,  int _iColsOut,  double *_pdblRcond)
564 {
565     int iReturn = 0;
566     int iIndex  = 0;
567     char cNorm  = 0;
568     int iExit   = 0;
569     int iLower = 0;
570     int iUpper = 0;
571
572     /*temporary variables*/
573     int iWorkMin    = 0;
574     int iInfo       = 0;
575     int iMax        = 0;
576     double dblRcond = 0;
577
578     double dblEps       = 0;
579     double RCONDthresh  = 0;
580     double dblAnorm     = 0;
581     double dblOne       = 1;
582
583     double *pAf     = NULL;
584     double *pXb     = NULL;
585     double *pDwork  = NULL;
586
587     double* dblTemp = NULL;
588     int iOne        = 1;
589     int iSize       = 0;
590
591     int *pRank  = NULL;
592     int *pIpiv  = NULL;
593     int *pJpvt  = NULL;
594     int *pIwork = NULL;
595
596     iWorkMin    = Max(4 * _iCols1, Max(Min(_iRows1, _iCols1) + 3 * _iCols1 + 1, 2 * Min(_iRows1, _iCols1) + _iCols2));
597
598     /* Array allocations*/
599     pAf         = (double*)malloc(sizeof(double) * _iRows1 * _iCols1);
600     pXb         = (double*)malloc(sizeof(double) * Max(_iRows1, _iCols1) * _iCols2);
601
602     pRank       = (int*)malloc(sizeof(int));
603     pIpiv       = (int*)malloc(sizeof(int) * _iCols1);
604     pJpvt       = (int*)malloc(sizeof(int) * _iCols1);
605     pIwork      = (int*)malloc(sizeof(int) * _iCols1);
606
607     pDwork      = (double*)malloc(sizeof(double) * iWorkMin);
608     dblEps      = nc_eps();
609     RCONDthresh = 10 * dblEps;
610
611     if (_iRows1 == _iCols1)
612     {
613         matrixIsTriangular(_pdblReal1, NULL, _iRows1, _iCols1, &iUpper, &iLower);
614         if (iUpper || iLower)
615         {
616             //if matrix is triangular
617             char cUpLoType[4] = {'N','U','L','U'};
618             char cUpLo = cUpLoType[iUpper + 2*iLower];
619             C2F(dtrcon)("1", &cUpLo, "N", &_iRows1, _pdblReal1, &_iRows1, &dblRcond, pDwork, pIwork, &iInfo);
620             if (dblRcond > RCONDthresh)
621             {
622                 // _pdblReal2 will be overwritten by dtrsm
623                 iSize = _iRows2 * _iCols2;
624                 dblTemp = (double*)malloc(iSize * sizeof(double));
625                 C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne);
626
627                 //solve triangular system
628                 C2F(dtrsm)("L", &cUpLo,"N", "N", &_iRows2, &_iCols2, &dblOne, _pdblReal1, &_iRows1, dblTemp, &_iRows2);
629                 C2F(dcopy)(&iSize, dblTemp, &iOne, _pdblRealOut, &iOne);
630                 iExit = 1;
631
632                 free(dblTemp);
633             }
634         }
635         else
636         {
637             iSize = _iRows1 * _iCols1;
638             C2F(dcopy)(&iSize, _pdblReal1, &iOne, pAf, &iOne);
639             C2F(dgetrf)(&_iCols1, &_iCols1, pAf, &_iCols1, pIpiv, &iInfo);
640             if (iInfo == 0)
641             {
642                 dblAnorm = C2F(dlange)("1", &_iRows1, &_iCols1, _pdblReal1, &_iRows1, pDwork);
643                 C2F(dgecon)("1", &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
644                 if (dblRcond > RCONDthresh)
645                 {
646                     // _pdblReal2 will be overwritten by dgetrs
647                     iSize = _iRows2 * _iCols2;
648                     dblTemp = (double*)malloc(iSize * sizeof(double));
649                     C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne);
650
651                     C2F(dgetrs)("N", &_iCols1, &_iCols2, pAf, &_iCols1, pIpiv, dblTemp, &_iCols1, &iInfo);
652                     C2F(dcopy)(&iSize, dblTemp, &iOne, _pdblRealOut, &iOne);
653                     iExit = 1;
654
655                     free(dblTemp);
656                 }
657             }
658         }
659
660         if (iExit == 0)
661         {
662             *_pdblRcond = dblRcond;
663             iReturn = -1;
664         }
665     }
666
667     if (iExit == 0)
668     {
669         dblRcond = RCONDthresh;
670         cNorm = 'F';
671         iMax = Max(_iRows1, _iCols1);
672         C2F(dlacpy)(&cNorm, &_iRows1, &_iCols2, _pdblReal2, &_iRows1, pXb, &iMax);
673         memset(pJpvt, 0x00, sizeof(int) * _iCols1);
674         // _pdblReal1 will be overwritten by dgelsy1
675         iSize = _iRows1 * _iCols1;
676         dblTemp = (double*)malloc(iSize * sizeof(double));
677         C2F(dcopy)(&iSize, _pdblReal1, &iOne, dblTemp, &iOne);
678         iInfo = 1;
679         C2F(dgelsy1)(&_iRows1, &_iCols1, &_iCols2, dblTemp, &_iRows1, pXb, &iMax,
680                      pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo);
681         free(dblTemp);
682
683         if (iInfo == 0)
684         {
685             if ( _iRows1 != _iCols1 && pRank[0] < Min(_iRows1, _iCols1))
686             {
687                 iReturn = -2;
688                 *_pdblRcond = pRank[0];
689             }
690
691             cNorm = 'F';
692             C2F(dlacpy)(&cNorm, &_iCols1, &_iCols2, pXb, &iMax, _pdblRealOut, &_iCols1);
693         }
694     }
695
696     free(pAf);
697     free(pXb);
698     free(pRank);
699     free(pIpiv);
700     free(pJpvt);
701     free(pIwork);
702     free(pDwork);
703     return 0;
704 }
705
706
707 /*Complex matrices left division*/
708 int iLeftDivisionOfComplexMatrix(
709     double *_pdblReal1,     double *_pdblImg1,      int _iRows1,    int _iCols1,
710     double *_pdblReal2,     double *_pdblImg2,      int _iRows2,    int _iCols2,
711     double *_pdblRealOut,   double *_pdblImgOut,    int _iRowsOut,  int _iColsOut,  double *_pdblRcond)
712 {
713     int iReturn = 0;
714     int iIndex  = 0;
715     char cNorm  = 0;
716     int iExit   = 0;
717     int iLower = 0;
718     int iUpper = 0;
719
720     /*temporary variables*/
721     int iWorkMin    = 0;
722     int iInfo       = 0;
723     int iMax        = 0;
724     double dblRcond = 0;
725
726     double dblEps       = 0;
727     double RCONDthresh  = 0;
728     double dblAnorm     = 0;
729     doublecomplex dblCplxOne = {1,0};
730
731     doublecomplex *pAf      = NULL;
732     doublecomplex *pXb      = NULL;
733     doublecomplex *pDwork   = NULL;
734     doublecomplex *poVar1   = NULL;
735     doublecomplex *poVar2   = NULL;
736     doublecomplex *poOut    = NULL;
737
738     double *pRwork  = NULL;
739
740     int iRank       = 0;
741     int *pIpiv      = NULL;
742     int *pJpvt      = NULL;
743
744     iWorkMin    = Max(2 * _iCols1, Min(_iRows1, _iCols1) + Max(2 * Min(_iRows1, _iCols1), Max(_iCols1, Min(_iRows1, _iCols1) + _iCols2)));
745
746     /* Array allocations*/
747     poVar1      = oGetDoubleComplexFromPointer(_pdblReal1, _pdblImg1, _iRows1 * _iCols1);
748     poVar2      = oGetDoubleComplexFromPointer(_pdblReal2, _pdblImg2, _iRows2 * _iCols2);
749
750     pIpiv       = (int*)malloc(sizeof(int) * _iCols1);
751     pJpvt       = (int*)malloc(sizeof(int) * _iCols1);
752     pRwork      = (double*)malloc(sizeof(double) * _iCols1 * 2);
753
754     cNorm       = '1';
755     pDwork      = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
756     dblEps      = nc_eps();
757     RCONDthresh = 10 * dblEps;
758
759     if (_iRows1 == _iCols1)
760     {
761         matrixIsTriangular(_pdblReal1, _pdblImg1, _iRows1, _iCols1, &iUpper, &iLower);
762         if (iUpper || iLower)
763         {
764             //if matrix is triangular
765             char cUpLoType[4] = {'N','U','L','U'};
766             char cUpLo = cUpLoType[iUpper + 2*iLower];
767             C2F(ztrcon)("1", &cUpLo, "N", &_iRows1, poVar1, &_iRows1, &dblRcond, pDwork, pRwork, &iInfo);
768             if (dblRcond > RCONDthresh)
769             {
770                 //solve triangular system
771                 C2F(ztrsm)("L", &cUpLo,"N", "N", &_iRows2, &_iCols2, &dblCplxOne, poVar1, &_iRows1, poVar2, &_iRows2);
772                 vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
773                 iExit = 1;
774             }
775             else
776             {
777                 //how to extract that ? Oo
778                 iReturn = -1;
779                 *_pdblRcond = dblRcond;
780             }
781         }
782         else
783         {
784             C2F(zgetrf)(&_iCols1, &_iCols1, poVar1, &_iCols1, pIpiv, &iInfo);
785             if (iInfo == 0)
786             {
787                 dblAnorm = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
788                 C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
789                 if (dblRcond > RCONDthresh)
790                 {
791                     C2F(zgetrs)("N", &_iCols1, &_iCols2, poVar1, &_iCols1, pIpiv, poVar2, &_iCols1, &iInfo);
792                     vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
793                     iExit = 1;
794                 }
795                 else
796                 {
797                     //how to extract that ? Oo
798                     iReturn = -1;
799                     *_pdblRcond = dblRcond;
800                 }
801             }
802         }
803     }
804
805     if (iExit == 0)
806     {
807         dblRcond = RCONDthresh;
808         iMax = Max(_iRows1, _iCols1);
809         memset(pJpvt, 0x00, sizeof(int) * _iCols1);
810         pXb = (doublecomplex*)malloc(sizeof(doublecomplex) * iMax * _iColsOut);
811         cNorm = 'F';
812         C2F(zlacpy)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)pXb, &iMax);
813         // pXb : in input pXb is of size rows1 x col2
814         //       in output pXp is of size col1 x col2
815         iInfo = 1;
816         C2F(zgelsy1)(&_iRows1, &_iCols1, &_iCols2, poVar1, &_iRows1, pXb, &iMax,
817                      pJpvt, &dblRcond, &iRank, pDwork, &iWorkMin, pRwork, &iInfo);
818
819         if (iInfo == 0)
820         {
821             // In the case where "pXb" has more rows that the output,
822             // the output values are the first lines of pXb
823             // and not the size of output first elements of pXb.
824             double* tmpRealPart = (double*)malloc(iMax * _iColsOut * sizeof(double));
825             double* tmpImagPart = (double*)malloc(iMax * _iColsOut * sizeof(double));
826             vGetPointerFromDoubleComplex(pXb, iMax * _iColsOut, tmpRealPart, tmpImagPart);
827
828             if ( _iRows1 != _iCols1 && iRank < Min(_iRows1, _iCols1))
829             {
830                 //how to extract that ? Oo
831                 iReturn = -2;
832                 *_pdblRcond = (double)iRank;
833             }
834
835             C2F(dlacpy)(&cNorm, &_iRowsOut, &_iColsOut, tmpRealPart, &iMax, _pdblRealOut, &_iRowsOut);
836             C2F(dlacpy)(&cNorm, &_iRowsOut, &_iColsOut, tmpImagPart, &iMax, _pdblImgOut, &_iRowsOut);
837
838             free(tmpRealPart);
839             free(tmpImagPart);
840         }
841         free(pXb);
842     }
843
844     vFreeDoubleComplexFromPointer(poVar1);
845     vFreeDoubleComplexFromPointer(poVar2);
846     free(pIpiv);
847     free(pJpvt);
848     free(pRwork);
849     free(pDwork);
850     return 0;
851 }
852
853 void matrixIsTriangular(double *_pdblReal, double *_pdblImg, int _iRows, int _iCols, int *bUpper, int *bLower)
854 {
855     double *pdblArray[2] = {_pdblReal, _pdblImg};
856     double *pdbl;
857
858     *bUpper = 1;
859     *bLower = 1;
860
861     for (int i=0; i<2 && pdblArray[i] != NULL; i++)
862     {
863         int j;
864         int iDim;
865         int iOne = 1;
866         //upper triangular ?
867         if (*bUpper)
868         {
869             pdbl = pdblArray[i] + 1;
870             iDim = _iRows-1;
871             for (j = 0; j < _iCols; j++, iDim--)
872             {
873                 // compute L1 norm with BLAS appears to be faster than testing
874                 // successive nullity of individual terms
875                 if (C2F(dasum)(&iDim, pdbl, &iOne) > 0.0) break;
876                 pdbl += _iRows+1;
877             }
878             *bUpper = (j==_iCols);
879         }
880         //lower triangular ?
881         if (*bLower)
882         {
883             pdbl = pdblArray[i] + _iRows;
884             for (j = 1; j < _iCols; j++)
885             {
886                 if (C2F(dasum)(&j, pdbl, &iOne) > 0.0) break;
887                 pdbl += _iRows;
888             }
889             *bLower = (j==_iCols);
890         }
891     }
892 }