linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / src / c / svd.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2009 - DIGITEO - Bernard HUGUENEY
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 #include <stdlib.h>
13 #include <string.h>
14 #include <stdio.h>
15 #include "doublecomplex.h"
16 #include "machine.h"
17 #include "MALLOC.h"
18 #include "core_math.h"
19 #include "svd.h"
20
21 /*
22  * Lapack functions performing the real work
23  */
24 extern void C2F(zgesvd)(char const jobU[1]/* 'A'|'S'|'O'|'N'*/, char const jobVT[1]/* 'A'|'S'|'O'|'N'*/, int const* iRows, int const* iCols
25                    , doublecomplex* a, int const* ldA, double* s, doublecomplex* u, int const* ldU, doublecomplex* vt, int const* ldVT
26                    , doublecomplex* work, int const* lWork, double* rWork, int* info);
27 extern void C2F(dgesvd)(char const jobU[1]/* 'A'|'S'|'O'|'N'*/, char const jobVT[1]/* 'A'|'S'|'O'|'N'*/, int const* iRows, int const* iCols
28                    , double* a, int const* ldA, double* s, double* u, int const* ldU, double* vt, int const* ldVT
29                    , double* work, int const* lWork, int* info);
30
31 /*
32  * Lapack utility functions
33  */
34 /* sorting (used to correct a Lapack bug */
35 extern void C2F(dlasrt)(char const id[1]/* 'I'ncreasing|'D'ecreasing order */, int const* length, double* array, int* info);
36 /* querying for epsilon, used for default treshold value */
37 extern double C2F(dlamch)(char const query[1], long );
38
39 /*
40  * functions used to allocate workspace for Lapack functions
41  * in: job char 
42  *
43  */
44 static double* allocDgesvdWorkspace(char job, int iRows, int iCols, int iColsToCompute, int* iAllocated);
45 static int iDgesvd(char job, double* pData, int iRows, int iCols, int iColsToCompute, double* pSV, double* pU, double* pVT, double* pWork, int iWorkSize);
46 static doublecomplex* allocZgesvdWorkspace(char job, int iRows, int iCols, int computeA, int* iAllocated);
47 static int iZgesvd(char job, doublecomplex* pData, int iRows, int iCols, int iColsToCompute, double* pSV, doublecomplex* pU, doublecomplex* pVT
48                    , doublecomplex* pWork, int iWorkSize, double* pRWork);
49
50 /* original fortran code says this is needed to workaround a bug in lapack, should assess necessity for current Lapack*/
51 static int iCorrectSV(double* pSV, int iMinRowsCols);
52 /* if lhs >1, S and V are computed from SV and VT */
53 static int iComputeSandVReal( double const* pSV, double const* pVT, int iRows, int iCols, int economy, double* pS, double* pV);
54 static int iComputeSandVComplex( double const* pSV, doublecomplex const* pVT, int iRows, int iCols, int economy, double* pS, doublecomplex* pV);
55
56
57 static doublecomplex* allocZgesvdWorkspace(char job, int iRows, int iCols, int iColsToCompute, int* iAllocated)
58 {
59   doublecomplex* res= NULL;
60   doublecomplex optimal;
61   int query= -1; /* query zgesvd for optimal worksize */
62   int info;
63   C2F(zgesvd)(&job, &job, &iRows, &iCols, NULL, &iRows, NULL, NULL, &iRows, NULL, &iColsToCompute, &optimal, &query, NULL, &info);
64   *iAllocated= (int) optimal.r;
65   res= MALLOC((*iAllocated) * sizeof(doublecomplex)); /* try to allocate optimal worksize */
66   if(!res) /* allocation failed, try with minimal worksize */
67     {
68       *iAllocated= 2 * Min(iRows, iCols) + Max(iRows, iCols);
69       res= MALLOC((*iAllocated) * sizeof(doublecomplex));
70       if(!res) /* failed again, report that we allocated 0 bytes */
71         {
72           *iAllocated= 0;
73         }
74     }
75   return res;
76 }
77 static int iZgesvd(char job, doublecomplex* pData, int iRows, int iCols, int iColsToCompute, double* pSV, doublecomplex* pU, doublecomplex* pVT
78                    , doublecomplex* pWork, int iWorkSize, double* pRWork)
79 {
80   int info;
81   C2F(zgesvd)(&job, &job, &iRows, &iCols, pData, &iRows, pSV, (pU ? pU : pData) , &iRows, (pVT ? pVT : pData), &iColsToCompute
82               , pWork, &iWorkSize, pRWork, &info);
83   return ( job == 'N' ) ? iCorrectSV(pSV, Min(iRows, iCols)) : info ;
84 }
85
86 static double* allocDgesvdWorkspace(char job, int iRows, int iCols, int iColsToCompute, int* iAllocated)
87 {
88   double* res=NULL;
89   double optimal;
90   int query=-1;
91   int info;
92   C2F(dgesvd)(&job, &job, &iRows, &iCols, NULL, &iRows, NULL, NULL, &iRows, NULL, &iColsToCompute, &optimal, &query, &info);
93   *iAllocated= (int) optimal;
94   res= MALLOC((*iAllocated) * sizeof(double));
95   if(!res)
96     {
97       *iAllocated= Max(3*Min(iRows, iCols) + Max(iRows, iCols), 5* Min(iRows, iCols)); /* Max & Min might have been more eficient as functions */
98       res= MALLOC((*iAllocated) * sizeof(double));
99       if(!res)
100         {
101           *iAllocated= 0;
102         }
103     }
104   return res;
105 }
106
107 static int iDgesvd(char job, double* pData, int iRows, int iCols, int iColsToCompute, double* pSV, double* pU, double* pVT
108                    , double* pWork, int iWorkSize)
109 {
110   int info;
111   C2F(dgesvd)(&job, &job, &iRows, &iCols, pData, &iRows, pSV, (pU ? pU : pData) , &iRows, (pVT ? pVT : pData), &iColsToCompute
112               , pWork, &iWorkSize, &info);
113   return ( job == 'N' ) ? iCorrectSV(pSV, Min(iRows, iCols)) : info ;
114 }
115
116 /*
117   lhs==1 ->   pSV && !pU && !pS && !pV
118   lhs==3|4 ->  !pSV && pU && pS && pV
119   lhs==4 <-> pRk != NULL
120  */
121
122 int iSvdM(double* pData, int iRows, int iCols, int complexArg, int economy, double tol, double* pSV, double* pU, double* pS, double* pV, double* pRk)
123 {
124   int ret=0, allocOK;
125   int worksize;
126
127   double* pWork= NULL;
128   double* pRWork= NULL;
129   double* pVT= NULL;
130
131   char const job= pU ? (economy ? 'S' :'A') : 'N';
132   int colsToCompute= (job=='S') ? Min(iRows, iCols) : iCols;
133   pSV= pSV ? pSV : (double*) MALLOC(Min(iRows, iCols) * sizeof(double));
134   if( (allocOK= ( (pWork= (complexArg 
135                           ? (double*)allocZgesvdWorkspace(job, iRows, iCols, colsToCompute, &worksize)
136                           : allocDgesvdWorkspace(job, iRows, iCols, colsToCompute, &worksize)))
137                  && (!complexArg 
138                      || (pRWork= (double*) MALLOC(5 * Min(iRows, iCols) * sizeof(double))))
139                  && (pVT= (pU 
140                            ? (double*)MALLOC(colsToCompute * iCols * (complexArg ? sizeof(doublecomplex) : sizeof(double)))
141                            : pData)))))
142     {
143       ret= complexArg 
144         ? iZgesvd(job,  (doublecomplex*)pData, iRows, iCols, colsToCompute, pSV,(doublecomplex*) pU, (doublecomplex*)pVT
145                   , (doublecomplex*)pWork, worksize, pRWork)
146         : iDgesvd(job, pData, iRows, iCols, colsToCompute, pSV, pU, pVT, pWork, worksize);
147       if( (ret==0) && pU)
148         {
149           /* openmp sections */
150           /* openmp section */
151           {
152             int unused= complexArg 
153               ? iComputeSandVComplex(pSV, (doublecomplex*)pVT, iRows, iCols, economy, pS, (doublecomplex*)pV)
154               : iComputeSandVReal(pSV, pVT, iRows, iCols, economy, pS, pV);
155           }
156           /* openmp section */
157           {
158             if(pRk) /* compute rank */
159               {
160                 int i, rk;
161                 tol = (tol == 0.) ? (Max(iRows, iCols) * C2F(dlamch)("e",1L) * pSV[0]) : tol; /* original Fortran code does the strict fp compare */
162                 rk= -1;
163                 for(i=0; i!= Min(iRows, iCols); ++i)
164                 {
165                   if(pSV[i] > tol)
166                     {
167                       rk= i;
168                     }
169                 }
170                 *pRk= (double)(rk+1);
171               }
172           }
173         }
174     }
175   FREE(pWork);
176   FREE(pRWork);
177   if(pU)
178     {
179       FREE(pVT);
180     }
181   return allocOK ? ret : -1;
182 }
183
184
185 static int iCorrectSV(double* pSV, int iMinRowsCols)
186 {
187   int i, info;
188   for(i=0; i != iMinRowsCols; ++i)
189     {
190       pSV[i]= Abs(pSV[i]); /* should be fabs in c99 */
191     }
192   C2F(dlasrt)("D", &iMinRowsCols, pSV, &info);
193   return info;
194 }
195
196 static int iComputeSandVComplex( double const* pSV, doublecomplex const* pVT, int iRows, int iCols, int economy, double* pS, doublecomplex* pV)
197 {
198   int i, j;
199   int minRowsCols= Min(iRows, iCols);
200   int economyCols= economy ? minRowsCols : iCols;
201   int economyRows= economy ? minRowsCols : iRows;
202   memset(pS, 0, economyRows * economyCols * sizeof(double));
203   for(i=0; i != minRowsCols; ++i)
204     {
205       pS[i+i*economyRows]= pSV[i];
206     }
207   /* multicore: with at least 2 caches, each half ( in !economy mode) should be done in it's on thread to avoid cache trashing */
208   for(j=0; j != economyCols; ++j)
209     {
210       for(i= 0; i != iCols; ++i)
211         {
212           pV[i + j*iCols].r= pVT[j + i*economyCols].r;
213           pV[i + j*iCols].i= -pVT[j + i*economyCols].i;
214           if(!economy)
215             {
216               pV[j + i*economyCols].r= pVT[i + j*iCols].r;
217               pV[j + i*economyCols].i= -pVT[i + j*iCols].i;
218             }
219         }
220     }
221   return 0;
222 }
223
224 static int iComputeSandVReal( double const* pSV, double const* pVT, int iRows, int iCols, int economy, double* pS, double* pV)
225 {
226   int i, j;
227   int minRowsCols= Min(iRows, iCols);
228   int economyCols= economy ? minRowsCols : iCols;
229   int economyRows= economy ? minRowsCols : iRows;
230
231   memset(pS, 0, economyRows * economyCols * sizeof(double));
232   for(i=0; i != minRowsCols; ++i)
233     {
234       pS[i+i*economyRows]= pSV[i];
235     }
236   for(j=0; j != economyCols; ++j)
237     {
238       for(i= 0; i != iCols; ++i)
239         {
240           pV[i + j*iCols]= pVT[j + i*economyCols];
241           if(!economy)
242             {
243               pV[j + i*iCols]= pVT[i + j*iCols];
244             }
245         }
246     }
247   return 0;
248 }