1 /*
2 Copyright (c) 2009-2014, Jack Poulson
3 All rights reserved.
4
5 This file is part of Elemental and is under the BSD 2-Clause License,
6 which can be found in the LICENSE file in the root directory, or at
7 http://opensource.org/licenses/BSD-2-Clause
8 */
9 #pragma once
10 #ifndef ELEM_SKELETON_HPP
11 #define ELEM_SKELETON_HPP
12
13 #include ELEM_ADJOINT_INC
14 #include ELEM_QR_BUSINGERGOLUB_INC
15 #include ELEM_PERMUTECOLS_INC
16 #include ELEM_PSEUDOINVERSE_INC
17
18 // NOTE: There are *many* algorithms for (pseudo-)skeleton/CUR decompositions,
19 // and, for now, we will simply implement one.
20
21 // TODO: More algorithms and more options (e.g., default tolerances).
22
23 // TODO: Implement randomized algorithms from Jiawei Chiu and Laurent Demanet's
24 // "Sublinear randomized algorithms for skeleton decompositions"?
25
26 namespace elem {
27
28 template<typename F>
29 inline void
Skeleton(const Matrix<F> & A,Matrix<Int> & permR,Matrix<Int> & permC,Matrix<F> & Z,Int maxSteps,Base<F> tol)30 Skeleton
31 ( const Matrix<F>& A,
32 Matrix<Int>& permR, Matrix<Int>& permC,
33 Matrix<F>& Z, Int maxSteps, Base<F> tol )
34 {
35 DEBUG_ONLY(CallStackEntry cse("Skeleton"))
36 // Find the row permutation
37 Matrix<F> B;
38 Adjoint( A, B );
39 const Int numSteps = qr::BusingerGolub( B, permR, maxSteps, tol );
40
41 // Form pinv(AR')=pinv(AR)'
42 Adjoint( A, B );
43 InversePermuteCols( B, permR );
44 B.Resize( B.Height(), numSteps );
45 Pseudoinverse( B );
46
47 // Form K := A pinv(AR)
48 Matrix<F> K;
49 Gemm( NORMAL, ADJOINT, F(1), A, B, K );
50
51 // Find the column permutation (force the same number of steps)
52 B = A;
53 qr::BusingerGolub( B, permC, numSteps );
54
55 // Form pinv(AC)
56 B = A;
57 InversePermuteCols( B, permC );
58 B.Resize( B.Height(), numSteps );
59 Pseudoinverse( B );
60
61 // Form Z := pinv(AC) K = pinv(AC) (A pinv(AR))
62 Gemm( NORMAL, NORMAL, F(1), B, K, Z );
63 }
64
65 template<typename F,Dist UPerm>
66 inline void
Skeleton(const DistMatrix<F> & A,DistMatrix<Int,UPerm,STAR> & permR,DistMatrix<Int,UPerm,STAR> & permC,DistMatrix<F> & Z,Int maxSteps,Base<F> tol)67 Skeleton
68 ( const DistMatrix<F>& A,
69 DistMatrix<Int,UPerm,STAR>& permR, DistMatrix<Int,UPerm,STAR>& permC,
70 DistMatrix<F>& Z, Int maxSteps, Base<F> tol )
71 {
72 DEBUG_ONLY(CallStackEntry cse("Skeleton"))
73 const Grid& g = A.Grid();
74
75 // Find the row permutation
76 DistMatrix<F> B(g);
77 Adjoint( A, B );
78 const Int numSteps = qr::BusingerGolub( B, permR, maxSteps, tol );
79
80 // Form pinv(AR')=pinv(AR)'
81 Adjoint( A, B );
82 InversePermuteCols( B, permR );
83 B.Resize( B.Height(), numSteps );
84 Pseudoinverse( B );
85
86 // Form K := A pinv(AR)
87 DistMatrix<F> K(g);
88 Gemm( NORMAL, ADJOINT, F(1), A, B, K );
89
90 // Find the column permutation (force the same number of steps)
91 B = A;
92 qr::BusingerGolub( B, permC, numSteps );
93
94 // Form pinv(AC)
95 B = A;
96 InversePermuteCols( B, permC );
97 B.Resize( B.Height(), numSteps );
98 Pseudoinverse( B );
99
100 // Form Z := pinv(AC) K = pinv(AC) (A pinv(AR))
101 Gemm( NORMAL, NORMAL, F(1), B, K, Z );
102 }
103
104 } // namespace elem
105
106 #endif // ifndef ELEM_SKELETON_HPP
107