1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup spop_strans
20 //! @{
21
22
23
24 template<typename eT>
25 inline
26 void
apply_noalias(SpMat<eT> & B,const SpMat<eT> & A)27 spop_strans::apply_noalias(SpMat<eT>& B, const SpMat<eT>& A)
28 {
29 arma_extra_debug_sigprint();
30
31 B.reserve(A.n_cols, A.n_rows, A.n_nonzero); // deliberately swapped
32
33 if(A.n_nonzero == 0) { return; }
34
35 // This follows the TRANSP algorithm described in
36 // 'Sparse Matrix Multiplication Package (SMMP)'
37 // (R.E. Bank and C.C. Douglas, 2001)
38
39 const uword m = A.n_rows;
40 const uword n = A.n_cols;
41
42 const eT* a = A.values;
43 eT* b = access::rwp(B.values);
44
45 const uword* ia = A.col_ptrs;
46 const uword* ja = A.row_indices;
47
48 uword* ib = access::rwp(B.col_ptrs);
49 uword* jb = access::rwp(B.row_indices);
50
51 // // ib is already zeroed, as B is freshly constructed
52 //
53 // for(uword i=0; i < (m+1); ++i)
54 // {
55 // ib[i] = 0;
56 // }
57
58 for(uword i=0; i < n; ++i)
59 {
60 for(uword j = ia[i]; j < ia[i+1]; ++j)
61 {
62 ib[ ja[j] + 1 ]++;
63 }
64 }
65
66 for(uword i=0; i < m; ++i)
67 {
68 ib[i+1] += ib[i];
69 }
70
71 for(uword i=0; i < n; ++i)
72 {
73 for(uword j = ia[i]; j < ia[i+1]; ++j)
74 {
75 const uword jj = ja[j];
76
77 const uword ib_jj = ib[jj];
78
79 jb[ib_jj] = i;
80
81 b[ib_jj] = a[j];
82
83 ib[jj]++;
84 }
85 }
86
87 for(uword i = m-1; i >= 1; --i)
88 {
89 ib[i] = ib[i-1];
90 }
91
92 ib[0] = 0;
93 }
94
95
96
97 template<typename T1>
98 inline
99 void
apply(SpMat<typename T1::elem_type> & out,const SpOp<T1,spop_strans> & in)100 spop_strans::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_strans>& in)
101 {
102 arma_extra_debug_sigprint();
103
104 typedef typename T1::elem_type eT;
105
106 const unwrap_spmat<T1> U(in.m);
107
108 if(U.is_alias(out))
109 {
110 SpMat<eT> tmp;
111
112 spop_strans::apply_noalias(tmp, U.M);
113
114 out.steal_mem(tmp);
115 }
116 else
117 {
118 spop_strans::apply_noalias(out, U.M);
119 }
120 }
121
122
123
124 //! for transpose of non-complex matrices, redirected from spop_htrans::apply()
125 template<typename T1>
126 inline
127 void
apply(SpMat<typename T1::elem_type> & out,const SpOp<T1,spop_htrans> & in)128 spop_strans::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_htrans>& in)
129 {
130 arma_extra_debug_sigprint();
131
132 typedef typename T1::elem_type eT;
133
134 const unwrap_spmat<T1> U(in.m);
135
136 if(U.is_alias(out))
137 {
138 SpMat<eT> tmp;
139
140 spop_strans::apply_noalias(tmp, U.M);
141
142 out.steal_mem(tmp);
143 }
144 else
145 {
146 spop_strans::apply_noalias(out, U.M);
147 }
148 }
149
150
151
152 //! @}
153