linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / src / c / rcond.c
1
2 /*
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) ????-2009 - INRIA
5 *
6 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which
8 * you should have received as part of this distribution.  The terms
9 * are also available at
10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13 #include "machine.h"
14 #include "doublecomplex.h"
15 #include "MALLOC.h"
16 #include "rcond.h"
17
18 /* Lapack routines doing the real work */
19 extern double C2F(zlange)(char const norm[1], int const* m, int const* n, doublecomplex const* a, int const* ldA, double* work, int* info);
20 extern double C2F(dlange)(char const norm[1], int const* m, int const* n, double const* a, int const* ldA, double* work, int* info);
21 extern void C2F(zgecon)(char const norm[1], int const* n, doublecomplex const* a, int const* ldA, double const* aNorm, double* rCond,
22     doublecomplex* work, double* rWork, int* info);
23 extern void C2F(dgecon)(char const norm[1], int const* n, double const* a, int const* ldA, double const* aNorm, double* rCond,
24     double* work, int* iWork, int* info);
25 extern void C2F(zgetrf)(int const* m, int const* n, doublecomplex* a, int const* ldA, int* iPiv, int* info);
26 extern void C2F(dgetrf)(int const* m, int const* n, double* a, int const* ldA, int* iPiv, int* info);
27
28 /* part of the API, cf. rcond.h */
29 int iRcondM(double* pData, int iCols, int complexArg, double* pRcond)
30 {
31     int info;
32     int* piPiv= (int*) MALLOC( iCols * sizeof(int));
33     void* pRIWork= complexArg ? MALLOC(2*iCols * sizeof(double)) : MALLOC(iCols * sizeof(int));
34     double* pWork= (double*)MALLOC(complexArg ? (2*iCols*sizeof(doublecomplex)) : (4*iCols*sizeof(double)));
35     if( piPiv && pRIWork && pWork)
36     {
37         double aNorm = 0;
38         if(complexArg)
39         {
40             aNorm = C2F(zlange)("1", &iCols, &iCols, (doublecomplex*)pData, &iCols, NULL, &info);
41         }
42         else
43         {
44             aNorm = C2F(dlange)("1", &iCols, &iCols, pData, &iCols, NULL, &info);
45         }
46
47         *pRcond= 0.;
48
49         if(complexArg)
50         {
51             C2F(zgetrf)(&iCols, &iCols, (doublecomplex*)pData, &iCols, piPiv, &info);
52             if(!info)
53             {
54                 C2F(zgecon)("1", &iCols, (doublecomplex*)pData, &iCols, &aNorm, pRcond, (doublecomplex*)pWork, (double*)pRIWork, &info);
55             }
56         }
57         else
58         {
59             C2F(dgetrf)(&iCols, &iCols, pData, &iCols, piPiv, &info);
60             if(!info)
61             {
62                 C2F(dgecon)("1", &iCols, pData, &iCols, &aNorm, pRcond, pWork, (int*)pRIWork, &info);
63             }
64         }
65     }
66     else
67     {
68         info= -1; /* MALLOC error */
69     }
70     FREE(piPiv);
71     FREE(pRIWork);
72     FREE(pWork);
73     return info;
74 }