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