1/* ffpack_rankprofiles.inl
2 * Copyright (C) 2015 FFLAS-FFACK group
3 *
4 * Written by Clement Pernet <Clement.Pernet@imag.fr>
5 *
6 * ========LICENCE========
7 * This file is part of the library FFLAS-FFPACK.
8 *
9 * FFLAS-FFPACK is free software: you can redistribute it and/or modify
10 * it under the terms of the  GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 * This library is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with this library; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22 * ========LICENCE========
23 *.
24 */
25
26#ifndef __FFLASFFPACK_ffpack_rank_profiles_INL
27#define __FFLASFFPACK_ffpack_rank_profiles_INL
28
29namespace FFPACK{
30    template <class Field>
31    inline size_t RowRankProfile (const Field& F, const size_t M, const size_t N,
32                                  typename Field::Element_ptr A, const size_t lda,
33                                  size_t* &rkprofile, const FFPACK_LU_TAG LuTag){
34        FFLAS::ParSeqHelper::Sequential seqH;
35        return FFPACK::RowRankProfile (F, M, N, A, lda, rkprofile, LuTag, seqH);
36    }
37
38    template <class Field>
39    inline size_t pRowRankProfile (const Field& F, const size_t M, const size_t N,
40                                  typename Field::Element_ptr A, const size_t lda,
41                                   size_t* &rkprofile, size_t numthreads, const FFPACK_LU_TAG LuTag){
42        size_t r;
43        PAR_BLOCK{
44            size_t nt = numthreads ? numthreads : NUM_THREADS;
45            FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::Threads> parH(nt);
46            r = FFPACK::RowRankProfile (F, M, N, A, lda, rkprofile, LuTag, parH);
47        }
48        return r;
49    }
50
51    template <class Field, class PSHelper>
52    inline size_t RowRankProfile (const Field& F, const size_t M, const size_t N,
53                                  typename Field::Element_ptr A, const size_t lda,
54                                  size_t* &rkprofile, const FFPACK_LU_TAG LuTag, PSHelper& psH){
55
56
57        size_t *P = FFLAS::fflas_new<size_t>((LuTag==FfpackSlabRecursive)?N:M);
58        size_t *Q = FFLAS::fflas_new<size_t>((LuTag==FfpackSlabRecursive)?M:N);
59        size_t R;
60
61        if (LuTag == FfpackSlabRecursive){
62            R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, M, N, A, lda, P, Q);
63            std::swap(P,Q);
64        } else
65            R = PLUQ (F, FFLAS::FflasNonUnit, M, N, A, lda, P, Q, psH);
66
67        rkprofile = FFLAS::fflas_new<size_t> (R);
68
69        RankProfileFromLU (P, M, R, rkprofile, LuTag);
70
71        FFLAS::fflas_delete (Q);
72        FFLAS::fflas_delete (P);
73        return R;
74    }
75
76    template <class Field>
77    inline size_t ColumnRankProfile (const Field& F, const size_t M, const size_t N,
78                                     typename Field::Element_ptr A, const size_t lda,
79                                     size_t* &rkprofile, const FFPACK_LU_TAG LuTag){
80        FFLAS::ParSeqHelper::Sequential seqH;
81        return FFPACK::ColumnRankProfile (F, M, N, A, lda, rkprofile, LuTag, seqH);
82    }
83
84    template <class Field>
85    inline size_t pColumnRankProfile (const Field& F, const size_t M, const size_t N,
86                                     typename Field::Element_ptr A, const size_t lda,
87                                      size_t* &rkprofile, size_t numthreads, const FFPACK_LU_TAG LuTag){
88        size_t r;
89        PAR_BLOCK{
90            size_t nt = numthreads ? numthreads : NUM_THREADS;
91            FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::Threads> parH(nt);
92            r = FFPACK::ColumnRankProfile (F, M, N, A, lda, rkprofile, LuTag, parH);
93        }
94        return r;
95    }
96
97    template <class Field, class PSHelper>
98    inline size_t ColumnRankProfile (const Field& F, const size_t M, const size_t N,
99                              typename Field::Element_ptr A, const size_t lda,
100                              size_t* &rkprofile, const FFPACK_LU_TAG LuTag, PSHelper& psH){
101
102        size_t *P = FFLAS::fflas_new<size_t>(M);
103        size_t *Q = FFLAS::fflas_new<size_t>(N);
104        size_t R;
105
106        if (LuTag == FfpackSlabRecursive){
107            R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasTrans, M, N, A, lda, P, Q);
108        } else
109            R = PLUQ (F, FFLAS::FflasNonUnit, M, N, A, lda, P, Q, psH);
110
111        rkprofile = FFLAS::fflas_new<size_t> (R);
112
113        RankProfileFromLU (Q, N, R, rkprofile, LuTag);
114
115        FFLAS::fflas_delete (P);
116        FFLAS::fflas_delete (Q);
117        return R;
118    }
119
120
121    inline void RankProfileFromLU (const size_t* Q, const size_t N, const size_t R,
122                                   size_t* rkprofile, const FFPACK_LU_TAG LuTag){
123
124        if (LuTag == FfpackSlabRecursive)
125            std::copy(Q, Q+R, rkprofile);
126        else {
127            size_t * RP = FFLAS::fflas_new<size_t>(N);
128            for (size_t i=0;i < N; ++i)
129                RP [i] = i;
130            for (size_t i=0; i<N; ++i)
131                if (Q[i] != i)
132                    std::swap (RP [i], RP [Q [i]]);
133
134            std::copy(RP, RP+R, rkprofile);
135            std::sort (rkprofile, rkprofile + R);
136            FFLAS::fflas_delete(RP);
137        }
138    }
139
140    inline size_t LeadingSubmatrixRankProfiles (const size_t M, const size_t N, const size_t R,
141                                                const size_t LSm, const size_t LSn,
142                                                const size_t* P, const size_t* Q,
143                                                size_t* RRP, size_t* CRP){
144        size_t LSr=0; // rank of the LSm x LSn leading submatrix
145
146        size_t* MathP = FFLAS::fflas_new<size_t>(M);
147        size_t* MathQ = FFLAS::fflas_new<size_t>(N);
148
149        LAPACKPerm2MathPerm (MathP, P, M);
150        LAPACKPerm2MathPerm (MathQ, Q, N);
151        for (size_t i = 0; i < R; i++)
152            if (MathP[i] < LSm && MathQ[i] < LSn){
153                RRP [LSr] = MathP[i];
154                CRP [LSr] = MathQ[i];
155                LSr++;
156            }
157        std::sort (RRP, RRP+LSr);
158        std::sort (CRP, CRP+LSr);
159        FFLAS::fflas_delete(MathP);
160        FFLAS::fflas_delete(MathQ);
161        return LSr;
162
163    }
164
165
166    template <class Field>
167    size_t RowRankProfileSubmatrixIndices (const Field& F,
168                                           const size_t M, const size_t N,
169                                           typename Field::Element_ptr A,
170                                           const size_t lda,
171                                           size_t*& rowindices,
172                                           size_t*& colindices,
173                                           size_t& R)
174    {
175        size_t *P = FFLAS::fflas_new<size_t>(N);
176        size_t *Q = FFLAS::fflas_new<size_t>(M);
177
178        R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, M, N, A, lda, P, Q);
179        rowindices = FFLAS::fflas_new<size_t>(M);
180        colindices = FFLAS::fflas_new<size_t>(N);
181        for (size_t i=0; i<R; ++i){
182            rowindices [i] = Q [i];
183        }
184        for (size_t i=0; i<N; ++i)
185            colindices [i] = i;
186        size_t tmp;
187        for (size_t i=0; i<R; ++i){
188            if (i != P[i]){
189                tmp = colindices[i];
190                colindices[i] = colindices[P[i]];
191                colindices[P[i]] = tmp;
192            }
193        }
194
195        FFLAS::fflas_delete( P);
196        FFLAS::fflas_delete( Q);
197
198        return R;
199    }
200
201    template <class Field>
202    size_t ColRankProfileSubmatrixIndices (const Field& F,
203                                           const size_t M, const size_t N,
204                                           typename Field::Element_ptr A,
205                                           const size_t lda,
206                                           size_t*& rowindices,
207                                           size_t*& colindices,
208                                           size_t& R)
209    {
210        size_t *P = FFLAS::fflas_new<size_t>(M);
211        size_t *Q = FFLAS::fflas_new<size_t>(N);
212
213        R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasTrans, M, N, A, lda, P, Q);
214        rowindices = FFLAS::fflas_new<size_t>(M);
215        colindices = FFLAS::fflas_new<size_t>(N);
216        for (size_t i=0; i<R; ++i)
217            colindices [i] = Q [i];
218
219        for (size_t i=0; i<N; ++i)
220            rowindices [i] = i;
221
222        size_t tmp;
223        for (size_t i=0; i<R; ++i){
224            if (i != P[i]){
225                tmp = rowindices[i];
226                rowindices[i] = rowindices[P[i]];
227                rowindices[P[i]] = tmp;
228            }
229        }
230        FFLAS::fflas_delete( P);
231        FFLAS::fflas_delete( Q);
232
233        return R;
234    }
235
236    template <class Field>
237    size_t RowRankProfileSubmatrix (const Field& F,
238                                    const size_t M, const size_t N,
239                                    typename Field::Element_ptr A,
240                                    const size_t lda,
241                                    typename Field::Element_ptr& X, size_t& R)
242    {
243
244        size_t * rowindices, * colindices;
245
246        typename Field::Element_ptr A2 = FFLAS::fflas_new (F, M, N) ;
247        FFLAS::fassign(F,M,N,A,lda,A2,N);
248
249        RowRankProfileSubmatrixIndices (F, M, N, A2, N, rowindices, colindices, R);
250
251        X = FFLAS::fflas_new (F, R, R);
252        for (size_t i=0; i<R; ++i)
253            for (size_t j=0; j<R; ++j)
254                F.assign (*(X + i*R + j), *(A + rowindices[i]*lda + colindices[j]));
255        FFLAS::fflas_delete (A2);
256        FFLAS::fflas_delete( rowindices);
257        FFLAS::fflas_delete( colindices);
258        return R;
259    }
260
261    template <class Field>
262    size_t ColRankProfileSubmatrix (const Field& F, const size_t M, const size_t N,
263                                    typename Field::Element_ptr A, const size_t lda,
264                                    typename Field::Element_ptr& X, size_t& R)
265    {
266
267        size_t * rowindices, * colindices;
268
269        typename Field::Element_ptr A2 = FFLAS::fflas_new (F, M, N);
270        FFLAS::fassign(F,M,N,A,lda,A2,N);
271
272        ColRankProfileSubmatrixIndices (F, M, N, A2, N, rowindices, colindices, R);
273
274        X = FFLAS::fflas_new (F, R, R);
275        for (size_t i=0; i<R; ++i)
276            for (size_t j=0; j<R; ++j)
277                F.assign (*(X + i*R + j), *(A + rowindices[i]*lda + colindices[j]));
278        FFLAS::fflas_delete (A2);
279        FFLAS::fflas_delete( colindices);
280        FFLAS::fflas_delete( rowindices);
281        return R;
282    }
283
284    template <class Field>
285    typename Field::Element_ptr
286    LQUPtoInverseOfFullRankMinor( const Field& F, const size_t rank,
287                                  typename Field::Element_ptr A_factors, const size_t lda,
288                                  const size_t* QtPointer,
289                                  typename Field::Element_ptr X, const size_t ldx)
290    {
291
292        // upper entries are okay, just need to move up bottom ones
293        const size_t* srcRow = QtPointer;
294        for (size_t row=0; row<rank; row++, srcRow++)
295            if (*srcRow != row) {
296                typename Field::Element_ptr oldRow = A_factors + (*srcRow) * lda;
297                typename Field::Element_ptr newRow = A_factors + row * lda;
298                for (size_t col=0; col<row; col++, oldRow++, newRow++)
299                    F.assign(*newRow, *oldRow);
300            }
301
302        // X <- (Qt.L.Q)^(-1)
303        //invL( F, rank, A_factors, lda, X, ldx);
304        ftrtri (F, FFLAS::FflasLower, FFLAS::FflasUnit, rank, A_factors, lda);
305        FFLAS::fassign(F,rank,rank,A_factors,lda,X,ldx);
306
307        // X = U^-1.X
308        ftrsm( F, FFLAS::FflasLeft, FFLAS::FflasUpper, FFLAS::FflasNoTrans,
309               FFLAS::FflasNonUnit, rank, rank, F.one, A_factors, lda, X, ldx);
310
311        return X;
312
313    }
314
315} // namespace FFPACK
316
317#endif // __FFLASFFPACK_ffpack_rank_profiles_INL
318/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
319// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
320