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