summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorCedric Delamarre <cedric.delamarre@scilab-enterprises.com>2013-07-02 12:13:02 +0200
committerCedric Delamarre <cedric.delamarre@scilab-enterprises.com>2013-07-02 12:15:02 +0200
commit0a585148ce38e9c420bf78f5a4ae7eae4bd09a00 (patch)
tree49712c8e9a7a0e68bd0fad7031759179b7194669
parent7407373dbe97d0b98f6336f1a5b664f686e68774 (diff)
downloadscilab-YaSp2012.zip
scilab-YaSp2012.tar.gz
matrix division updated.YaSp2012
After master modification : https://codereview.scilab.org/#/c/11610 https://codereview.scilab.org/#/c/10705 https://codereview.scilab.org/#/c/10753 test_run("linear_algebra","bug_9196",["no_check_error_output" ]); Change-Id: I9a7bc445ab14a8e0359d1b8d0ad332e7fc0d0a10
-rw-r--r--scilab/modules/operations/src/c/matrix_division.c426
1 files changed, 218 insertions, 208 deletions
diff --git a/scilab/modules/operations/src/c/matrix_division.c b/scilab/modules/operations/src/c/matrix_division.c
index 2bf7029..0331369 100644
--- a/scilab/modules/operations/src/c/matrix_division.c
+++ b/scilab/modules/operations/src/c/matrix_division.c
@@ -16,15 +16,15 @@
16 16
17/*iRightDivisionComplexMatrixByComplexMatrix*/ 17/*iRightDivisionComplexMatrixByComplexMatrix*/
18int iRightDivisionComplexMatrixByComplexMatrix( 18int iRightDivisionComplexMatrixByComplexMatrix(
19 double *_pdblReal1, double *_pdblImg1, int _iInc1, 19 double *_pdblReal1, double *_pdblImg1, int _iInc1,
20 double *_pdblReal2, double *_pdblImg2, int _iInc2, 20 double *_pdblReal2, double *_pdblImg2, int _iInc2,
21 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) 21 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize)
22{ 22{
23 int iErr = 0; 23 int iErr = 0;
24 int iIndex = 0; //Main loop index 24 int iIndex = 0; //Main loop index
25 int iIndex1 = 0; //Loop index on left operand 25 int iIndex1 = 0; //Loop index on left operand
26 int iIndex2 = 0; //Loop index on right operand 26 int iIndex2 = 0; //Loop index on right operand
27 int iIndexOut = 0; //Lopp index on result matrix 27 int iIndexOut = 0; //Lopp index on result matrix
28 28
29 if (_iInc2 == 0) 29 if (_iInc2 == 0)
30 { 30 {
@@ -37,9 +37,9 @@ int iRightDivisionComplexMatrixByComplexMatrix(
37 for (iIndex = 0 ; iIndex < _iSize ; iIndex++) 37 for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
38 { 38 {
39 iErr = iRightDivisionComplexByComplex(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); 39 iErr = iRightDivisionComplexByComplex(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
40 iIndexOut += _iIncOut; 40 iIndexOut += _iIncOut;
41 iIndex1 += _iInc1; 41 iIndex1 += _iInc1;
42 iIndex2 += _iInc2; 42 iIndex2 += _iInc2;
43 } 43 }
44 44
45 return iErr; 45 return iErr;
@@ -58,19 +58,19 @@ int iRightDivisionComplexByComplex(
58 { 58 {
59 //got NaN + i NaN 59 //got NaN + i NaN
60 iErr = 10; 60 iErr = 10;
61 *_pdblRealOut = _dblImg2 / _dblReal2; 61 *_pdblRealOut = _dblImg2 / _dblReal2;
62 *_pdblImgOut = *_pdblRealOut; 62 *_pdblImgOut = *_pdblRealOut;
63 } 63 }
64 else 64 else
65 { 65 {
66 *_pdblRealOut = _dblReal1 / _dblReal2; 66 *_pdblRealOut = _dblReal1 / _dblReal2;
67 *_pdblImgOut = _dblImg1 / _dblReal2; 67 *_pdblImgOut = _dblImg1 / _dblReal2;
68 } 68 }
69 } 69 }
70 else if (_dblReal2 == 0) 70 else if (_dblReal2 == 0)
71 { 71 {
72 *_pdblRealOut = _dblImg1 / _dblImg2; 72 *_pdblRealOut = _dblImg1 / _dblImg2;
73 *_pdblImgOut = (-_dblReal1) / _dblImg2; 73 *_pdblImgOut = (-_dblReal1) / _dblImg2;
74 } 74 }
75 else 75 else
76 { 76 {
@@ -78,17 +78,17 @@ int iRightDivisionComplexByComplex(
78 78
79 if (dabss(_dblReal2) >= dabss(_dblImg2)) 79 if (dabss(_dblReal2) >= dabss(_dblImg2))
80 { 80 {
81 double dblRatio2 = _dblImg2 / _dblReal2; 81 double dblRatio2 = _dblImg2 / _dblReal2;
82 double dblSum = _dblReal2 + dblRatio2 * _dblImg2; 82 double dblSum = _dblReal2 + dblRatio2 * _dblImg2;
83 *_pdblRealOut = (_dblReal1 + _dblImg1 * dblRatio2) / dblSum; 83 *_pdblRealOut = (_dblReal1 + _dblImg1 * dblRatio2) / dblSum;
84 *_pdblImgOut = (_dblImg1 - _dblReal1 * dblRatio2) / dblSum; 84 *_pdblImgOut = (_dblImg1 - _dblReal1 * dblRatio2) / dblSum;
85 } 85 }
86 else 86 else
87 { 87 {
88 double dblRatio2 = _dblReal2 / _dblImg2; 88 double dblRatio2 = _dblReal2 / _dblImg2;
89 double dblSum = _dblImg2 + dblRatio2 * _dblReal2; 89 double dblSum = _dblImg2 + dblRatio2 * _dblReal2;
90 *_pdblRealOut = (_dblReal1 * dblRatio2 + _dblImg1) / dblSum; 90 *_pdblRealOut = (_dblReal1 * dblRatio2 + _dblImg1) / dblSum;
91 *_pdblImgOut = (_dblImg1 * dblRatio2 - _dblReal1) / dblSum; 91 *_pdblImgOut = (_dblImg1 * dblRatio2 - _dblReal1) / dblSum;
92 } 92 }
93 } 93 }
94 return iErr; 94 return iErr;
@@ -96,15 +96,15 @@ int iRightDivisionComplexByComplex(
96 96
97/*iRightDivisionRealMatrixByComplexMatrix*/ 97/*iRightDivisionRealMatrixByComplexMatrix*/
98int iRightDivisionRealMatrixByComplexMatrix( 98int iRightDivisionRealMatrixByComplexMatrix(
99 double *_pdblReal1, int _iInc1, 99 double *_pdblReal1, int _iInc1,
100 double *_pdblReal2, double *_pdblImg2, int _iInc2, 100 double *_pdblReal2, double *_pdblImg2, int _iInc2,
101 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) 101 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize)
102{ 102{
103 int iErr = 0; 103 int iErr = 0;
104 int iIndex = 0; //Main loop index 104 int iIndex = 0; //Main loop index
105 int iIndex1 = 0; //Loop index on left operand 105 int iIndex1 = 0; //Loop index on left operand
106 int iIndex2 = 0; //Loop index on right operand 106 int iIndex2 = 0; //Loop index on right operand
107 int iIndexOut = 0; //Lopp index on result matrix 107 int iIndexOut = 0; //Lopp index on result matrix
108 108
109 if (_iInc2 == 0) 109 if (_iInc2 == 0)
110 { 110 {
@@ -117,9 +117,9 @@ int iRightDivisionRealMatrixByComplexMatrix(
117 for (iIndex = 0 ; iIndex < _iSize ; iIndex++) 117 for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
118 { 118 {
119 iErr = iRightDivisionRealByComplex(_pdblReal1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); 119 iErr = iRightDivisionRealByComplex(_pdblReal1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
120 iIndexOut += _iIncOut; 120 iIndexOut += _iIncOut;
121 iIndex1 += _iInc1; 121 iIndex1 += _iInc1;
122 iIndex2 += _iInc2; 122 iIndex2 += _iInc2;
123 } 123 }
124 124
125 return iErr; 125 return iErr;
@@ -134,13 +134,13 @@ int iRightDivisionRealByComplex(
134 int iErr = 0; 134 int iErr = 0;
135 if (_dblImg2 == 0) 135 if (_dblImg2 == 0)
136 { 136 {
137 *_pdblRealOut = _dblReal1 / _dblReal2; 137 *_pdblRealOut = _dblReal1 / _dblReal2;
138 *_pdblImgOut = 0; 138 *_pdblImgOut = 0;
139 } 139 }
140 else if (_dblReal2 == 0) 140 else if (_dblReal2 == 0)
141 { 141 {
142 *_pdblRealOut = 0; 142 *_pdblRealOut = 0;
143 *_pdblImgOut = -_dblReal1 / _dblImg2; 143 *_pdblImgOut = -_dblReal1 / _dblImg2;
144 } 144 }
145 else 145 else
146 { 146 {
@@ -149,17 +149,17 @@ int iRightDivisionRealByComplex(
149 if (dblAbsSum == 0) 149 if (dblAbsSum == 0)
150 { 150 {
151 iErr = 10; 151 iErr = 10;
152 *_pdblRealOut = _dblReal1 / dblAbsSum; 152 *_pdblRealOut = _dblReal1 / dblAbsSum;
153 *_pdblImgOut = 0; 153 *_pdblImgOut = 0;
154 } 154 }
155 else 155 else
156 { 156 {
157 double dblReal1Sum = _dblReal1 / dblAbsSum; 157 double dblReal1Sum = _dblReal1 / dblAbsSum;
158 double dblReal2Sum = _dblReal2 / dblAbsSum; 158 double dblReal2Sum = _dblReal2 / dblAbsSum;
159 double dblImg2Sum = _dblImg2 / dblAbsSum; 159 double dblImg2Sum = _dblImg2 / dblAbsSum;
160 double dblSum = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2); 160 double dblSum = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2);
161 *_pdblRealOut = (dblReal1Sum * dblReal2Sum) / dblSum; 161 *_pdblRealOut = (dblReal1Sum * dblReal2Sum) / dblSum;
162 *_pdblImgOut = (-dblReal1Sum * dblImg2Sum) / dblSum; 162 *_pdblImgOut = (-dblReal1Sum * dblImg2Sum) / dblSum;
163 } 163 }
164 } 164 }
165 return iErr; 165 return iErr;
@@ -167,21 +167,21 @@ int iRightDivisionRealByComplex(
167 167
168/*iRightDivisionComplexMatrixByRealMatrix*/ 168/*iRightDivisionComplexMatrixByRealMatrix*/
169int iRightDivisionComplexMatrixByRealMatrix( 169int iRightDivisionComplexMatrixByRealMatrix(
170 double *_pdblReal1, double *_pdblImg1, int _iInc1, 170 double *_pdblReal1, double *_pdblImg1, int _iInc1,
171 double *_pdblReal2, int _iInc2, 171 double *_pdblReal2, int _iInc2,
172 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) 172 double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize)
173{ 173{
174 int iErr = 0; 174 int iErr = 0;
175 int iIndex = 0; //Main loop index 175 int iIndex = 0; //Main loop index
176 int iIndex1 = 0; //Loop index on left operand 176 int iIndex1 = 0; //Loop index on left operand
177 int iIndex2 = 0; //Loop index on right operand 177 int iIndex2 = 0; //Loop index on right operand
178 int iIndexOut = 0; //Lopp index on result matrix 178 int iIndexOut = 0; //Lopp index on result matrix
179 for (iIndex = 0 ; iIndex < _iSize ; iIndex++) 179 for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
180 { 180 {
181 iErr = iRightDivisionComplexByReal(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); 181 iErr = iRightDivisionComplexByReal(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]);
182 iIndexOut += _iIncOut; 182 iIndexOut += _iIncOut;
183 iIndex1 += _iInc1; 183 iIndex1 += _iInc1;
184 iIndex2 += _iInc2; 184 iIndex2 += _iInc2;
185 } 185 }
186 return iErr; 186 return iErr;
187} 187}
@@ -207,8 +207,8 @@ int iRightDivisionComplexByReal(
207 } 207 }
208 } 208 }
209 209
210 *_pdblRealOut = _dblReal1 / _dblReal2; 210 *_pdblRealOut = _dblReal1 / _dblReal2;
211 *_pdblImgOut = _dblImg1 / _dblReal2; 211 *_pdblImgOut = _dblImg1 / _dblReal2;
212 212
213 return iErr; 213 return iErr;
214} 214}
@@ -219,11 +219,11 @@ int iRightDivisionRealMatrixByRealMatrix(
219 double *_pdblReal2, int _iInc2, 219 double *_pdblReal2, int _iInc2,
220 double *_pdblRealOut, int _iIncOut, int _iSize) 220 double *_pdblRealOut, int _iIncOut, int _iSize)
221{ 221{
222 int iIndex = 0; //Main loop index 222 int iIndex = 0; //Main loop index
223 int iIndex1 = 0; //Loop index on left operand 223 int iIndex1 = 0; //Loop index on left operand
224 int iIndex2 = 0; //Loop index on right operand 224 int iIndex2 = 0; //Loop index on right operand
225 int iIndexOut = 0; //Lopp index on result matrix 225 int iIndexOut = 0; //Lopp index on result matrix
226 int iErr = 0; 226 int iErr = 0;
227 227
228 for (iIndex = 0 ; iIndex < _iSize ; iIndex++) 228 for (iIndex = 0 ; iIndex < _iSize ; iIndex++)
229 { 229 {
@@ -241,62 +241,64 @@ int iRightDivisionRealMatrixByRealMatrix(
241 } 241 }
242 242
243 _pdblRealOut[iIndexOut] = _pdblReal1[iIndex1] / _pdblReal2[iIndex2]; 243 _pdblRealOut[iIndexOut] = _pdblReal1[iIndex1] / _pdblReal2[iIndex2];
244 iIndexOut += _iIncOut; 244 iIndexOut += _iIncOut;
245 iIndex1 += _iInc1; 245 iIndex1 += _iInc1;
246 iIndex2 += _iInc2; 246 iIndex2 += _iInc2;
247 } 247 }
248 return iErr; 248 return iErr;
249} 249}
250 250
251int iRightDivisionOfRealMatrix( 251int iRightDivisionOfRealMatrix(
252 double *_pdblReal1, int _iRows1, int _iCols1, 252 double *_pdblReal1, int _iRows1, int _iCols1,
253 double *_pdblReal2, int _iRows2, int _iCols2, 253 double *_pdblReal2, int _iRows2, int _iCols2,
254 double *_pdblRealOut, int _iRowsOut, int _iColsOut, double* _pdblRcond) 254 double *_pdblRealOut, int _iRowsOut, int _iColsOut, double* _pdblRcond)
255{ 255{
256 int iReturn = 0; 256 int iReturn = 0;
257 257
258 int iIndex = 0; 258 int iIndex = 0;
259 char cNorm = 0; 259 char cNorm = 0;
260 int iExit = 0; 260 int iExit = 0;
261 261
262 /*temporary variables*/ 262 /*temporary variables*/
263 int iWorkMin = 0; 263 int iWorkMin = 0;
264 int iInfo = 0; 264 int iInfo = 0;
265 int iMax = 0; 265 int iMax = 0;
266 double dblRcond = 0; 266 double dblRcond = 0;
267 267
268 double dblEps = 0; 268 double dblEps = 0;
269 double dblAnorm = 0; 269 double RCONDthresh = 0;
270 double dblAnorm = 0;
270 271
271 double *pAf = NULL; 272 double *pAf = NULL;
272 double *pAt = NULL; 273 double *pAt = NULL;
273 double *pBt = NULL; 274 double *pBt = NULL;
274 double *pDwork = NULL; 275 double *pDwork = NULL;
275 276
276 int *pRank = NULL; 277 int *pRank = NULL;
277 int *pIpiv = NULL; 278 int *pIpiv = NULL;
278 int *pJpvt = NULL; 279 int *pJpvt = NULL;
279 int *pIwork = NULL; 280 int *pIwork = NULL;
280 281
281 iWorkMin = Max(4 * _iCols2, Max(Min(_iRows2, _iCols2) + 3 * _iRows2 + 1, 2 * Min(_iRows2, _iCols2) + _iRows1)); 282 iWorkMin = Max(4 * _iCols2, Max(Min(_iRows2, _iCols2) + 3 * _iRows2 + 1, 2 * Min(_iRows2, _iCols2) + _iRows1));
282 283
283 /* Array allocations*/ 284 /* Array allocations*/
284 pAf = (double*)malloc(sizeof(double) * _iCols2 * _iRows2); 285 pAf = (double*)malloc(sizeof(double) * _iCols2 * _iRows2);
285 pAt = (double*)malloc(sizeof(double) * _iCols2 * _iRows2); 286 pAt = (double*)malloc(sizeof(double) * _iCols2 * _iRows2);
286 pBt = (double*)malloc(sizeof(double) * Max(_iRows2, _iCols2) * _iRows1); 287 pBt = (double*)malloc(sizeof(double) * Max(_iRows2, _iCols2) * _iRows1);
287 288
288 pRank = (int*)malloc(sizeof(int)); 289 pRank = (int*)malloc(sizeof(int));
289 pIpiv = (int*)malloc(sizeof(int) * _iCols2); 290 pIpiv = (int*)malloc(sizeof(int) * _iCols2);
290 pJpvt = (int*)malloc(sizeof(int) * _iRows2); 291 pJpvt = (int*)malloc(sizeof(int) * _iRows2);
291 pIwork = (int*)malloc(sizeof(int) * _iCols2); 292 pIwork = (int*)malloc(sizeof(int) * _iCols2);
292 293
293 294
294 //C'est du grand nawak ca, on reserve toute la stack ! Oo 295 //C'est du grand nawak ca, on reserve toute la stack ! Oo
295 296
296 cNorm = '1'; 297 cNorm = '1';
297 pDwork = (double*)malloc(sizeof(double) * iWorkMin); 298 pDwork = (double*)malloc(sizeof(double) * iWorkMin);
298 dblEps = F2C(dlamch)("eps", 1L); 299 dblEps = F2C(dlamch)("eps", 1L);
299 dblAnorm = C2F(dlange)(&cNorm, &_iRows2, &_iCols2, _pdblReal2, &_iRows2, pDwork); 300 RCONDthresh = 10 * dblEps;
301 dblAnorm = C2F(dlange)(&cNorm, &_iRows2, &_iCols2, _pdblReal2, &_iRows2, pDwork);
300 302
301 //tranpose A and B 303 //tranpose A and B
302 304
@@ -310,23 +312,23 @@ int iRightDivisionOfRealMatrix(
310 { 312 {
311 ij = i + j * Max(_iRows2, _iCols2); 313 ij = i + j * Max(_iRows2, _iCols2);
312 ji = j + i * _iRows1; 314 ji = j + i * _iRows1;
313 pBt[ij] = _pdblReal1[ji]; 315 pBt[ij] = _pdblReal1[ji];
314 }//for(j = 0 ; j < _iRows1 ; j++) 316 }//for(j = 0 ; j < _iRows1 ; j++)
315 }//for(i = 0 ; i < _iCols2 ; i++) 317 }//for(i = 0 ; i < _iCols2 ; i++)
316 }//bloc esthetique 318 }//bloc esthetique
317 319
318 if (_iRows2 == _iCols2) 320 if (_iRows2 == _iCols2)
319 { 321 {
320 cNorm = 'F'; 322 cNorm = 'F';
321 C2F(dlacpy)(&cNorm, &_iCols2, &_iCols2, pAt, &_iCols2, pAf, &_iCols2); 323 C2F(dlacpy)(&cNorm, &_iCols2, &_iCols2, pAt, &_iCols2, pAf, &_iCols2);
322 C2F(dgetrf)(&_iCols2, &_iCols2, pAf, &_iCols2, pIpiv, &iInfo); 324 C2F(dgetrf)(&_iCols2, &_iCols2, pAf, &_iCols2, pIpiv, &iInfo);
323 if (iInfo == 0) 325 if (iInfo == 0)
324 { 326 {
325 cNorm = '1'; 327 cNorm = '1';
326 C2F(dgecon)(&cNorm, &_iCols2, pAf, &_iCols2, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo); 328 C2F(dgecon)(&cNorm, &_iCols2, pAf, &_iCols2, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
327 if (dblRcond > dsqrts(dblEps)) 329 if (dblRcond > RCONDthresh)
328 { 330 {
329 cNorm = 'N'; 331 cNorm = 'N';
330 C2F(dgetrs)(&cNorm, &_iCols2, &_iRows1, pAf, &_iCols2, pIpiv, pBt, &_iCols2, &iInfo); 332 C2F(dgetrs)(&cNorm, &_iCols2, &_iRows1, pAf, &_iCols2, pIpiv, pBt, &_iCols2, &iInfo);
331 vTransposeRealMatrix(pBt, _iCols2, _iRows1, _pdblRealOut); 333 vTransposeRealMatrix(pBt, _iCols2, _iRows1, _pdblRealOut);
332 iExit = 1; 334 iExit = 1;
@@ -343,10 +345,11 @@ int iRightDivisionOfRealMatrix(
343 345
344 if (iExit == 0) 346 if (iExit == 0)
345 { 347 {
346 dblRcond = dsqrts(dblEps); 348 dblRcond = RCONDthresh;
347 cNorm = 'F'; 349 cNorm = 'F';
348 iMax = Max(_iRows2, _iCols2); 350 iMax = Max(_iRows2, _iCols2);
349 memset(pJpvt, 0x00, sizeof(int) * _iRows2); 351 memset(pJpvt, 0x00, sizeof(int) * _iRows2);
352 iInfo = 1;
350 C2F(dgelsy1)(&_iCols2, &_iRows2, &_iRows1, pAt, &_iCols2, pBt, &iMax, 353 C2F(dgelsy1)(&_iCols2, &_iRows2, &_iRows1, pAt, &_iCols2, pBt, &iMax,
351 pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo); 354 pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo);
352 355
@@ -359,7 +362,7 @@ int iRightDivisionOfRealMatrix(
359 *_pdblRcond = pRank[0]; 362 *_pdblRcond = pRank[0];
360 } 363 }
361 364
362 // TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2); 365 // TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2);
363 366
364 //Mega caca de la mort qui tue des ours a mains nues 367 //Mega caca de la mort qui tue des ours a mains nues
365 //mais je ne sais pas comment le rendre "beau" :( 368 //mais je ne sais pas comment le rendre "beau" :(
@@ -371,7 +374,7 @@ int iRightDivisionOfRealMatrix(
371 { 374 {
372 ij = i + j * _iRows1; 375 ij = i + j * _iRows1;
373 ji = j + i * Max(_iRows2, _iCols2); 376 ji = j + i * Max(_iRows2, _iCols2);
374 _pdblRealOut[ij] = pBt[ji]; 377 _pdblRealOut[ij] = pBt[ji];
375 }//for(i = 0 ; i < _iRows2 ; i++) 378 }//for(i = 0 ; i < _iRows2 ; i++)
376 }//for(j = 0 ; j < _iRows1 ; j++) 379 }//for(j = 0 ; j < _iRows1 ; j++)
377 }//bloc esthetique 380 }//bloc esthetique
@@ -390,59 +393,61 @@ int iRightDivisionOfRealMatrix(
390} 393}
391 394
392 395
393int iRightDivisionOfComplexMatrix( 396int iRightDivisionOfComplexMatrix(
394 double *_pdblReal1, double *_pdblImg1, int _iRows1, int _iCols1, 397 double *_pdblReal1, double *_pdblImg1, int _iRows1, int _iCols1,
395 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2, 398 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2,
396 double *_pdblRealOut, double *_pdblImgOut, int _iRowsOut, int _iColsOut, double *_pdblRcond) 399 double *_pdblRealOut, double *_pdblImgOut, int _iRowsOut, int _iColsOut, double *_pdblRcond)
397{ 400{
398 int iReturn = 0; 401 int iReturn = 0;
399 int iIndex1 = 0; 402 int iIndex1 = 0;
400 int iIndex2 = 0; 403 int iIndex2 = 0;
401 char cNorm = 0; 404 char cNorm = 0;
402 int iExit = 0; 405 int iExit = 0;
403 406
404 /*temporary variables*/ 407 /*temporary variables*/
405 int iWorkMin = 0; 408 int iWorkMin = 0;
406 int iInfo = 0; 409 int iInfo = 0;
407 int iMax = 0; 410 int iMax = 0;
408 double dblRcond = 0; 411 double dblRcond = 0;
409 412
410 double dblEps = 0; 413 double dblEps = 0;
411 double dblAnorm = 0; 414 double RCONDthresh = 0;
415 double dblAnorm = 0;
412 416
413 doublecomplex *poVar1 = NULL; 417 doublecomplex *poVar1 = NULL;
414 doublecomplex *poVar2 = NULL; 418 doublecomplex *poVar2 = NULL;
415 doublecomplex *poOut = NULL; 419 doublecomplex *poOut = NULL;
416 doublecomplex *poAf = NULL; 420 doublecomplex *poAf = NULL;
417 doublecomplex *poAt = NULL; 421 doublecomplex *poAt = NULL;
418 doublecomplex *poBt = NULL; 422 doublecomplex *poBt = NULL;
419 doublecomplex *poDwork = NULL; 423 doublecomplex *poDwork = NULL;
420 424
421 int *pRank = NULL; 425 int *pRank = NULL;
422 int *pIpiv = NULL; 426 int *pIpiv = NULL;
423 int *pJpvt = NULL; 427 int *pJpvt = NULL;
424 double *pRwork = NULL; 428 double *pRwork = NULL;
425 429
426 iWorkMin = Max(2 * _iCols2, Min(_iRows2, _iCols2) + Max(2 * Min(_iRows2, _iCols2), Max(_iRows2 + 1, Min(_iRows2, _iCols2) + _iRows1))); 430 iWorkMin = Max(2 * _iCols2, Min(_iRows2, _iCols2) + Max(2 * Min(_iRows2, _iCols2), Max(_iRows2 + 1, Min(_iRows2, _iCols2) + _iRows1)));
427 431
428 /* Array allocations*/ 432 /* Array allocations*/
429 poVar1 = oGetDoubleComplexFromPointer(_pdblReal1, _pdblImg1, _iRows1 * _iCols1); 433 poVar1 = oGetDoubleComplexFromPointer(_pdblReal1, _pdblImg1, _iRows1 * _iCols1);
430 poVar2 = oGetDoubleComplexFromPointer(_pdblReal2, _pdblImg2, _iRows2 * _iCols2); 434 poVar2 = oGetDoubleComplexFromPointer(_pdblReal2, _pdblImg2, _iRows2 * _iCols2);
431 poOut = oGetDoubleComplexFromPointer(_pdblRealOut, _pdblImgOut, _iRowsOut * _iColsOut); 435 poOut = oGetDoubleComplexFromPointer(_pdblRealOut, _pdblImgOut, _iRowsOut * _iColsOut);
432 436
433 poAf = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2); 437 poAf = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2);
434 poAt = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2); 438 poAt = (doublecomplex*)malloc(sizeof(doublecomplex) * _iRows2 * _iCols2);
435 poBt = (doublecomplex*)malloc(sizeof(doublecomplex) * Max(_iRows2, _iCols2) * _iRows1); 439 poBt = (doublecomplex*)malloc(sizeof(doublecomplex) * Max(_iRows2, _iCols2) * _iRows1);
436 poDwork = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin); 440 poDwork = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
437 441
438 pRank = (int*)malloc(sizeof(int)); 442 pRank = (int*)malloc(sizeof(int));
439 pIpiv = (int*)malloc(sizeof(int) * _iCols2); 443 pIpiv = (int*)malloc(sizeof(int) * _iCols2);
440 pJpvt = (int*)malloc(sizeof(int) * _iRows2); 444 pJpvt = (int*)malloc(sizeof(int) * _iRows2);
441 pRwork = (double*)malloc(sizeof(double) * 2 * _iRows2); 445 pRwork = (double*)malloc(sizeof(double) * 2 * _iRows2);
442 446
443 dblEps = F2C(dlamch)("eps", 1L); 447 dblEps = F2C(dlamch)("eps", 1L);
444 cNorm = '1'; 448 RCONDthresh = 10 * dblEps;
445 dblAnorm = C2F(zlange)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)poDwork); 449 cNorm = '1';
450 dblAnorm = C2F(zlange)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)poDwork);
446 451
447 //tranpose A and B 452 //tranpose A and B
448 453
@@ -456,9 +461,9 @@ int iRightDivisionOfComplexMatrix(
456 { 461 {
457 ij = i + j * Max(_iRows2, _iCols2); 462 ij = i + j * Max(_iRows2, _iCols2);
458 ji = j + i * _iRows1; 463 ji = j + i * _iRows1;
459 poBt[ij].r = poVar1[ji].r; 464 poBt[ij].r = poVar1[ji].r;
460 //Conjugate 465 //Conjugate
461 poBt[ij].i = -poVar1[ji].i; 466 poBt[ij].i = -poVar1[ji].i;
462 }//for(j = 0 ; j < _iRows1 ; j++) 467 }//for(j = 0 ; j < _iRows1 ; j++)
463 }//for(i = 0 ; i < _iCols2 ; i++) 468 }//for(i = 0 ; i < _iCols2 ; i++)
464 }//bloc esthetique 469 }//bloc esthetique
@@ -466,16 +471,16 @@ int iRightDivisionOfComplexMatrix(
466 471
467 if (_iRows2 == _iCols2) 472 if (_iRows2 == _iCols2)
468 { 473 {
469 cNorm = 'F'; 474 cNorm = 'F';
470 C2F(zlacpy)(&cNorm, &_iCols2, &_iCols2, (double*)poAt, &_iCols2, (double*)poAf, &_iCols2); 475 C2F(zlacpy)(&cNorm, &_iCols2, &_iCols2, (double*)poAt, &_iCols2, (double*)poAf, &_iCols2);
471 C2F(zgetrf)(&_iCols2, &_iCols2, (double*)poAf, &_iCols2, pIpiv, &iInfo); 476 C2F(zgetrf)(&_iCols2, &_iCols2, (double*)poAf, &_iCols2, pIpiv, &iInfo);
472 if (iInfo == 0) 477 if (iInfo == 0)
473 { 478 {
474 cNorm = '1'; 479 cNorm = '1';
475 C2F(zgecon)(&cNorm, &_iCols2, poAf, &_iCols2, &dblAnorm, &dblRcond, poDwork, pRwork, &iInfo); 480 C2F(zgecon)(&cNorm, &_iCols2, poAf, &_iCols2, &dblAnorm, &dblRcond, poDwork, pRwork, &iInfo);
476 if (dblRcond > dsqrts(dblEps)) 481 if (dblRcond > RCONDthresh)
477 { 482 {
478 cNorm = 'N'; 483 cNorm = 'N';
479 C2F(zgetrs)(&cNorm, &_iCols2, &_iRows1, (double*)poAf, &_iCols2, pIpiv, (double*)poBt, &_iCols2, &iInfo); 484 C2F(zgetrs)(&cNorm, &_iCols2, &_iRows1, (double*)poAf, &_iCols2, pIpiv, (double*)poBt, &_iCols2, &iInfo);
480 vTransposeDoubleComplexMatrix(poBt, _iCols2, _iRows1, poOut, 1); 485 vTransposeDoubleComplexMatrix(poBt, _iCols2, _iRows1, poOut, 1);
481 vGetPointerFromDoubleComplex(poOut, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut); 486 vGetPointerFromDoubleComplex(poOut, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
@@ -493,10 +498,11 @@ int iRightDivisionOfComplexMatrix(
493 498
494 if (iExit == 0) 499 if (iExit == 0)
495 { 500 {
496 dblRcond = dsqrts(dblEps); 501 dblRcond = RCONDthresh;
497 cNorm = 'F'; 502 cNorm = 'F';
498 iMax = Max(_iRows2, _iCols2); 503 iMax = Max(_iRows2, _iCols2);
499 memset(pJpvt, 0x00, sizeof(int) * _iRows2); 504 memset(pJpvt, 0x00, sizeof(int) * _iRows2);
505 iInfo = 1;
500 C2F(zgelsy1)(&_iCols2, &_iRows2, &_iRows1, poAt, &_iCols2, poBt, &iMax, 506 C2F(zgelsy1)(&_iCols2, &_iRows2, &_iRows1, poAt, &_iCols2, poBt, &iMax,
501 pJpvt, &dblRcond, pRank, poDwork, &iWorkMin, pRwork, &iInfo); 507 pJpvt, &dblRcond, pRank, poDwork, &iWorkMin, pRwork, &iInfo);
502 508
@@ -509,7 +515,7 @@ int iRightDivisionOfComplexMatrix(
509 *_pdblRcond = pRank[0]; 515 *_pdblRcond = pRank[0];
510 } 516 }
511 517
512 // TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2); 518 // TransposeRealMatrix(pBt, _iRows1, _iRows2, _pdblRealOut, Max(_iRows1,_iCols1), _iRows2);
513 519
514 //Mega caca de la mort qui tue des ours a mains nues 520 //Mega caca de la mort qui tue des ours a mains nues
515 //mais je ne sais pas comment le rendre "beau" :( 521 //mais je ne sais pas comment le rendre "beau" :(
@@ -521,9 +527,9 @@ int iRightDivisionOfComplexMatrix(
521 { 527 {
522 ij = i + j * _iRows1; 528 ij = i + j * _iRows1;
523 ji = j + i * Max(_iRows2, _iCols2); 529 ji = j + i * Max(_iRows2, _iCols2);
524 _pdblRealOut[ij] = poBt[ji].r; 530 _pdblRealOut[ij] = poBt[ji].r;
525 //Conjugate 531 //Conjugate
526 _pdblImgOut[ij] = -poBt[ji].i; 532 _pdblImgOut[ij] = -poBt[ji].i;
527 }//for(i = 0 ; i < _iRows2 ; i++) 533 }//for(i = 0 ; i < _iRows2 ; i++)
528 }//for(j = 0 ; j < _iRows1 ; j++) 534 }//for(j = 0 ; j < _iRows1 ; j++)
529 }//bloc esthetique 535 }//bloc esthetique
@@ -547,84 +553,84 @@ int iRightDivisionOfComplexMatrix(
547} 553}
548 554
549/*Matrix left division*/ 555/*Matrix left division*/
550int iLeftDivisionOfRealMatrix( 556int iLeftDivisionOfRealMatrix(
551 double *_pdblReal1, int _iRows1, int _iCols1, 557 double *_pdblReal1, int _iRows1, int _iCols1,
552 double *_pdblReal2, int _iRows2, int _iCols2, 558 double *_pdblReal2, int _iRows2, int _iCols2,
553 double *_pdblRealOut, int _iRowsOut, int _iColsOut, double *_pdblRcond) 559 double *_pdblRealOut, int _iRowsOut, int _iColsOut, double *_pdblRcond)
554{ 560{
555 int iReturn = 0; 561 int iReturn = 0;
556 562 int iIndex = 0;
557 int iIndex = 0; 563 char cNorm = 0;
558 char cNorm = 0; 564 int iExit = 0;
559 int iExit = 0;
560 565
561 /*temporary variables*/ 566 /*temporary variables*/
562 int iWorkMin = 0; 567 int iWorkMin = 0;
563 int iInfo = 0; 568 int iInfo = 0;
564 int iMax = 0; 569 int iMax = 0;
565 double dblRcond = 0; 570 double dblRcond = 0;
566 571
567 double dblEps = 0; 572 double dblEps = 0;
568 double dblAnorm = 0; 573 double RCONDthresh = 0;
574 double dblAnorm = 0;
569 575
570 double *pAf = NULL; 576 double *pAf = NULL;
571 double *pXb = NULL; 577 double *pXb = NULL;
572 double *pDwork = NULL; 578 double *pDwork = NULL;
573 579
574 double* dblTemp = NULL; 580 double* dblTemp = NULL;
575 int iOne = 1; 581 int iOne = 1;
576 int iSize = 0; 582 int iSize = 0;
577 583
578 int *pRank = NULL; 584 int *pRank = NULL;
579 int *pIpiv = NULL; 585 int *pIpiv = NULL;
580 int *pJpvt = NULL; 586 int *pJpvt = NULL;
581 int *pIwork = NULL; 587 int *pIwork = NULL;
582 588
583 iWorkMin = Max(4 * _iCols1, Max(Min(_iRows1, _iCols1) + 3 * _iRows1 + 1, 2 * Min(_iRows1, _iCols1) + _iCols2)); 589 iWorkMin = Max(4 * _iCols1, Max(Min(_iRows1, _iCols1) + 3 * _iCols1 + 1, 2 * Min(_iRows1, _iCols1) + _iCols2));
584 590
585 /* Array allocations*/ 591 /* Array allocations*/
586 pAf = (double*)malloc(sizeof(double) * _iRows1 * _iCols1); 592 pAf = (double*)malloc(sizeof(double) * _iRows1 * _iCols1);
587 pXb = (double*)malloc(sizeof(double) * Max(_iRows1, _iCols1) * _iCols1); 593 pXb = (double*)malloc(sizeof(double) * Max(_iRows1, _iCols1) * _iCols2);
588 594
589 pRank = (int*)malloc(sizeof(int)); 595 pRank = (int*)malloc(sizeof(int));
590 pIpiv = (int*)malloc(sizeof(int) * _iCols1); 596 pIpiv = (int*)malloc(sizeof(int) * _iCols1);
591 pJpvt = (int*)malloc(sizeof(int) * _iCols1); 597 pJpvt = (int*)malloc(sizeof(int) * _iCols1);
592 pIwork = (int*)malloc(sizeof(int) * _iCols1); 598 pIwork = (int*)malloc(sizeof(int) * _iCols1);
593 599
594 cNorm = '1'; 600 cNorm = '1';
595 pDwork = (double*)malloc(sizeof(double) * iWorkMin); 601 pDwork = (double*)malloc(sizeof(double) * iWorkMin);
596 dblEps = F2C(dlamch)("eps", 1L); 602 dblEps = C2F(dlamch)("eps", 1L);
597 dblAnorm = C2F(dlange)(&cNorm, &_iRows1, &_iCols1, _pdblReal1, &_iRows1, pDwork); 603 RCONDthresh = 10 * dblEps;
604 dblAnorm = C2F(dlange)(&cNorm, &_iRows1, &_iCols1, _pdblReal1, &_iRows1, pDwork);
598 605
599 if (_iRows1 == _iCols1) 606 if (_iRows1 == _iCols1)
600 { 607 {
601 cNorm = 'F'; 608 cNorm = 'F';
602 C2F(dlacpy)(&cNorm, &_iCols1, &_iCols1, _pdblReal1, &_iCols1, pAf, &_iCols1); 609 C2F(dlacpy)(&cNorm, &_iCols1, &_iCols1, _pdblReal1, &_iCols1, pAf, &_iCols1);
603 C2F(dgetrf)(&_iCols1, &_iCols1, pAf, &_iCols1, pIpiv, &iInfo); 610 C2F(dgetrf)(&_iCols1, &_iCols1, pAf, &_iCols1, pIpiv, &iInfo);
604 if (iInfo == 0) 611 if (iInfo == 0)
605 { 612 {
606 cNorm = '1'; 613 cNorm = '1';
607 C2F(dgecon)(&cNorm, &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo); 614 C2F(dgecon)(&cNorm, &_iCols1, pAf, &_iCols1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
608 if (dblRcond > dsqrts(dblEps)) 615 if (dblRcond > RCONDthresh)
609 { 616 {
610 // _pdblReal2 will be overwrite by dgetrs 617 // _pdblReal2 will be overwrite by dgetrs
611 iSize = _iRows2 * _iCols2; 618 iSize = _iRows2 * _iCols2;
612 dblTemp = (double*)malloc(iSize * sizeof(double)); 619 dblTemp = (double*)malloc(iSize * sizeof(double));
613 C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne); 620 C2F(dcopy)(&iSize, _pdblReal2, &iOne, dblTemp, &iOne);
614 621
615 cNorm = 'N'; 622 cNorm = 'N';
616 C2F(dgetrs)(&cNorm, &_iCols1, &_iCols2, pAf, &_iCols1, pIpiv, dblTemp, &_iCols1, &iInfo); 623 C2F(dgetrs)(&cNorm, &_iCols1, &_iCols2, pAf, &_iCols1, pIpiv, dblTemp, &_iCols1, &iInfo);
617 cNorm = 'F'; 624 cNorm = 'F';
618 C2F(dlacpy)(&cNorm, &_iCols1, &_iCols2, dblTemp, &_iCols1, _pdblRealOut, &_iCols1); 625 C2F(dlacpy)(&cNorm, &_iCols1, &_iCols2, dblTemp, &_iCols1, _pdblRealOut, &_iCols1);
619 iExit = 1; 626 iExit = 1;
620 627
621 free(dblTemp); 628 free(dblTemp);
622 } 629 }
623 } 630 }
631
624 if (iExit == 0) 632 if (iExit == 0)
625 { 633 {
626 //how to extract that ? Oo
627 //how to extract that ? Oo
628 *_pdblRcond = dblRcond; 634 *_pdblRcond = dblRcond;
629 iReturn = -1; 635 iReturn = -1;
630 } 636 }
@@ -632,7 +638,7 @@ int iLeftDivisionOfRealMatrix(
632 638
633 if (iExit == 0) 639 if (iExit == 0)
634 { 640 {
635 dblRcond = dsqrts(dblEps); 641 dblRcond = RCONDthresh;
636 cNorm = 'F'; 642 cNorm = 'F';
637 iMax = Max(_iRows1, _iCols1); 643 iMax = Max(_iRows1, _iCols1);
638 C2F(dlacpy)(&cNorm, &_iRows1, &_iCols2, _pdblReal2, &_iRows1, pXb, &iMax); 644 C2F(dlacpy)(&cNorm, &_iRows1, &_iCols2, _pdblReal2, &_iRows1, pXb, &iMax);
@@ -641,14 +647,15 @@ int iLeftDivisionOfRealMatrix(
641 iSize = _iRows1 * _iCols1; 647 iSize = _iRows1 * _iCols1;
642 dblTemp = (double*)malloc(iSize * sizeof(double)); 648 dblTemp = (double*)malloc(iSize * sizeof(double));
643 C2F(dcopy)(&iSize, _pdblReal1, &iOne, dblTemp, &iOne); 649 C2F(dcopy)(&iSize, _pdblReal1, &iOne, dblTemp, &iOne);
644 C2F(dgelsy1)( &_iRows1, &_iCols1, &_iCols2, dblTemp, &_iRows1, pXb, &iMax, 650 iInfo = 1;
645 pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo); 651 C2F(dgelsy1)(&_iRows1, &_iCols1, &_iCols2, dblTemp, &_iRows1, pXb, &iMax,
652 pJpvt, &dblRcond, &pRank[0], pDwork, &iWorkMin, &iInfo);
646 free(dblTemp); 653 free(dblTemp);
654
647 if (iInfo == 0) 655 if (iInfo == 0)
648 { 656 {
649 if ( _iRows1 != _iCols1 && pRank[0] < Min(_iRows1, _iCols1)) 657 if ( _iRows1 != _iCols1 && pRank[0] < Min(_iRows1, _iCols1))
650 { 658 {
651 //how to extract that ? Oo
652 iReturn = -2; 659 iReturn = -2;
653 *_pdblRcond = pRank[0]; 660 *_pdblRcond = pRank[0];
654 } 661 }
@@ -686,8 +693,9 @@ int iLeftDivisionOfComplexMatrix(
686 int iMax = 0; 693 int iMax = 0;
687 double dblRcond = 0; 694 double dblRcond = 0;
688 695
689 double dblEps = 0; 696 double dblEps = 0;
690 double dblAnorm = 0; 697 double RCONDthresh = 0;
698 double dblAnorm = 0;
691 699
692 doublecomplex *pAf = NULL; 700 doublecomplex *pAf = NULL;
693 doublecomplex *pXb = NULL; 701 doublecomplex *pXb = NULL;
@@ -715,6 +723,7 @@ int iLeftDivisionOfComplexMatrix(
715 cNorm = '1'; 723 cNorm = '1';
716 pDwork = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin); 724 pDwork = (doublecomplex*)malloc(sizeof(doublecomplex) * iWorkMin);
717 dblEps = F2C(dlamch)("eps", 1L); 725 dblEps = F2C(dlamch)("eps", 1L);
726 RCONDthresh = 10 * dblEps;
718 dblAnorm = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork); 727 dblAnorm = C2F(zlange)(&cNorm, &_iRows1, &_iCols1, (double*)poVar1, &_iRows1, (double*)pDwork);
719 728
720 if (_iRows1 == _iCols1) 729 if (_iRows1 == _iCols1)
@@ -723,9 +732,9 @@ int iLeftDivisionOfComplexMatrix(
723 if (iInfo == 0) 732 if (iInfo == 0)
724 { 733 {
725 C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo); 734 C2F(zgecon)(&cNorm, &_iCols1, poVar1, &_iCols1, &dblAnorm, &dblRcond, pDwork, pRwork, &iInfo);
726 if (dblRcond > dsqrts(dblEps)) 735 if (dblRcond > RCONDthresh)
727 { 736 {
728 cNorm = 'N'; 737 cNorm = 'N';
729 C2F(zgetrs)(&cNorm, &_iCols1, &_iCols2, (double*)poVar1, &_iCols1, pIpiv, (double*)poVar2, &_iCols1, &iInfo); 738 C2F(zgetrs)(&cNorm, &_iCols1, &_iCols2, (double*)poVar1, &_iCols1, pIpiv, (double*)poVar2, &_iCols1, &iInfo);
730 vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut); 739 vGetPointerFromDoubleComplex(poVar2, _iRowsOut * _iColsOut, _pdblRealOut, _pdblImgOut);
731 iExit = 1; 740 iExit = 1;
@@ -741,7 +750,7 @@ int iLeftDivisionOfComplexMatrix(
741 750
742 if (iExit == 0) 751 if (iExit == 0)
743 { 752 {
744 dblRcond = dsqrts(dblEps); 753 dblRcond = RCONDthresh;
745 iMax = Max(_iRows1, _iCols1); 754 iMax = Max(_iRows1, _iCols1);
746 memset(pJpvt, 0x00, sizeof(int) * _iCols1); 755 memset(pJpvt, 0x00, sizeof(int) * _iCols1);
747 pXb = (doublecomplex*)malloc(sizeof(doublecomplex) * iMax * _iColsOut); 756 pXb = (doublecomplex*)malloc(sizeof(doublecomplex) * iMax * _iColsOut);
@@ -749,6 +758,7 @@ int iLeftDivisionOfComplexMatrix(
749 C2F(zlacpy)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)pXb, &iMax); 758 C2F(zlacpy)(&cNorm, &_iRows2, &_iCols2, (double*)poVar2, &_iRows2, (double*)pXb, &iMax);
750 // pXb : in input pXb is of size rows1 x col2 759 // pXb : in input pXb is of size rows1 x col2
751 // in output pXp is of size col1 x col2 760 // in output pXp is of size col1 x col2
761 iInfo = 1;
752 C2F(zgelsy1)(&_iRows1, &_iCols1, &_iCols2, poVar1, &_iRows1, pXb, &iMax, 762 C2F(zgelsy1)(&_iRows1, &_iCols1, &_iCols2, poVar1, &_iRows1, pXb, &iMax,
753 pJpvt, &dblRcond, &iRank, pDwork, &iWorkMin, pRwork, &iInfo); 763 pJpvt, &dblRcond, &iRank, pDwork, &iWorkMin, pRwork, &iInfo);
754 764