1 #ifndef CXXBLAS_INTERFACE_AUX_H
2 #define CXXBLAS_INTERFACE_AUX_H 1
3
4 #include <cxxblas/cxxblas.cxx>
5
6 using cxxblas::StorageOrder;
7 using cxxblas::ColMajor;
8 using cxxblas::RowMajor;
9 using cxxblas::StorageUpLo;
10 using cxxblas::Upper;
11 using cxxblas::Lower;
12
13
14 template <typename IndexType, typename MA>
15 void
16 switchFullStorageOrder(StorageOrder order, IndexType m, IndexType n,
17 const MA *&A, IndexType &ldA)
18 {
19 assert(order==ColMajor);
20
21 IndexType _ldA = n+1;
22 MA *_A = new MA[_ldA*m];
23
24 for (IndexType i=0; i<m; ++i) {
25 for (IndexType j=0; j<n; ++j) {
26 _A[i*_ldA+j] = A[j*ldA+i];
27 }
28 }
29
30 A = _A;
31 ldA = _ldA;
32 }
33
34 template <typename IndexType, typename MA>
35 void
36 switchFullStorageOrder(StorageOrder orderA, IndexType m, IndexType n,
37 const MA *A, IndexType &ldA,
38 MA *&_A, IndexType &_ldA)
39 {
40 if (orderA==RowMajor) {
41 for (IndexType i=0; i<m; ++i) {
42 for (IndexType j=0; j<n; ++j) {
43 _A[j*_ldA+i] = A[i*ldA+j];
44 }
45 }
46 } else {
47 for (IndexType i=0; i<m; ++i) {
48 for (IndexType j=0; j<n; ++j) {
49 _A[i*_ldA+j] = A[j*ldA+i];
50 }
51 }
52 }
53 }
54
55 template <typename IndexType, typename MA>
56 void
57 switchBandStorageOrder(StorageOrder order, IndexType m, IndexType n,
58 IndexType kl, IndexType ku,
59 const MA *&A, IndexType &ldA)
60 {
61 assert(order==ColMajor);
62
63 IndexType _ldA = kl+ku+1;
64 MA *_A = new MA[(kl+ku+1)*m];
65
66 for (IndexType i=0; i<m; ++i) {
67 for (IndexType j=0; j<n; ++j) {
68 if ((j-i<=ku) && (i-j<=kl)) {
69 _A[(kl+j-i)+_ldA*i] = A[(ku+i-j)+ldA*j];
70 }
71 }
72 }
73
74 A = _A;
75 ldA = _ldA;
76 }
77
78 template <typename IndexType, typename MA>
79 void
80 switchBandStorageOrder(StorageOrder order, StorageUpLo upLo,
81 IndexType n, IndexType k,
82 const MA *&A, IndexType &ldA)
83 {
84 IndexType kl = (upLo==Lower) ? k : 0;
85 IndexType ku = (upLo==Upper) ? k : 0;
86 switchBandStorageOrder(order, n, n, kl, ku, A, ldA);
87 }
88
89 template <typename IndexType, typename MA>
90 void
91 switchPackedStorageOrder(StorageOrder order, StorageUpLo &upLo,
92 IndexType n, const MA *&A)
93 {
94 assert(order==ColMajor);
95
96 MA *_A = new MA[n*(n+1)/2];
97 if (upLo==Upper) {
98 for (IndexType i=0; i<n; ++i) {
99 for (IndexType j=i; j<n; ++j) {
100 _A[j+i*(2*n-i-1)/2] = A[i+j*(j+1)/2];
101 }
102 }
103 } else {
104 for (IndexType i=0; i<n; ++i) {
105 for (IndexType j=0; j<=i; ++j) {
106 _A[j+i*(i+1)/2] = A[i+j*(2*n-j-1)/2];
107 }
108 }
109 }
110 A = _A;
111 }
112
113 template <typename IndexType, typename MA>
114 void
115 switchPackedStorageOrder(StorageOrder orderA, StorageUpLo &upLo,
116 IndexType n, const MA *A, MA *&_A)
117 {
118 if (upLo==Upper) {
119 if (orderA==ColMajor) {
120 for (IndexType i=0; i<n; ++i) {
121 for (IndexType j=i; j<n; ++j) {
122 _A[j+i*(2*n-i-1)/2] = A[i+j*(j+1)/2];
123 }
124 }
125 } else {
126 for (IndexType i=0; i<n; ++i) {
127 for (IndexType j=i; j<n; ++j) {
128 _A[i+j*(j+1)/2] = A[j+i*(2*n-i-1)/2];
129 }
130 }
131 }
132 } else {
133 if (orderA==ColMajor) {
134 for (IndexType i=0; i<n; ++i) {
135 for (IndexType j=0; j<=i; ++j) {
136 _A[j+i*(i+1)/2] = A[i+j*(2*n-j-1)/2];
137 }
138 }
139 } else {
140 for (IndexType i=0; i<n; ++i) {
141 for (IndexType j=0; j<=i; ++j) {
142 _A[i+j*(2*n-j-1)/2] = A[j+i*(i+1)/2];
143 }
144 }
145 }
146 }
147 }
148
149 template <typename IndexType, typename MA>
150 void
151 allocateFullStorage(StorageOrder order, IndexType m, IndexType n,
152 MA *&A, IndexType &ldA)
153 {
154 assert(order==RowMajor);
155
156 ldA = n;
157 A = new MA[ldA*m];
158 }
159
160 template <typename IndexType, typename MA>
161 void
162 allocatePackedStorage(IndexType n, MA *&A)
163 {
164 A = new MA[n*(n+1)/2];
165 }
166
167 template <typename MA>
168 void
169 releaseStorage(MA *A)
170 {
171 delete [] A;
172 }
173
174 #endif // CXXBLAS_INTERFACE_AUX_H
2 #define CXXBLAS_INTERFACE_AUX_H 1
3
4 #include <cxxblas/cxxblas.cxx>
5
6 using cxxblas::StorageOrder;
7 using cxxblas::ColMajor;
8 using cxxblas::RowMajor;
9 using cxxblas::StorageUpLo;
10 using cxxblas::Upper;
11 using cxxblas::Lower;
12
13
14 template <typename IndexType, typename MA>
15 void
16 switchFullStorageOrder(StorageOrder order, IndexType m, IndexType n,
17 const MA *&A, IndexType &ldA)
18 {
19 assert(order==ColMajor);
20
21 IndexType _ldA = n+1;
22 MA *_A = new MA[_ldA*m];
23
24 for (IndexType i=0; i<m; ++i) {
25 for (IndexType j=0; j<n; ++j) {
26 _A[i*_ldA+j] = A[j*ldA+i];
27 }
28 }
29
30 A = _A;
31 ldA = _ldA;
32 }
33
34 template <typename IndexType, typename MA>
35 void
36 switchFullStorageOrder(StorageOrder orderA, IndexType m, IndexType n,
37 const MA *A, IndexType &ldA,
38 MA *&_A, IndexType &_ldA)
39 {
40 if (orderA==RowMajor) {
41 for (IndexType i=0; i<m; ++i) {
42 for (IndexType j=0; j<n; ++j) {
43 _A[j*_ldA+i] = A[i*ldA+j];
44 }
45 }
46 } else {
47 for (IndexType i=0; i<m; ++i) {
48 for (IndexType j=0; j<n; ++j) {
49 _A[i*_ldA+j] = A[j*ldA+i];
50 }
51 }
52 }
53 }
54
55 template <typename IndexType, typename MA>
56 void
57 switchBandStorageOrder(StorageOrder order, IndexType m, IndexType n,
58 IndexType kl, IndexType ku,
59 const MA *&A, IndexType &ldA)
60 {
61 assert(order==ColMajor);
62
63 IndexType _ldA = kl+ku+1;
64 MA *_A = new MA[(kl+ku+1)*m];
65
66 for (IndexType i=0; i<m; ++i) {
67 for (IndexType j=0; j<n; ++j) {
68 if ((j-i<=ku) && (i-j<=kl)) {
69 _A[(kl+j-i)+_ldA*i] = A[(ku+i-j)+ldA*j];
70 }
71 }
72 }
73
74 A = _A;
75 ldA = _ldA;
76 }
77
78 template <typename IndexType, typename MA>
79 void
80 switchBandStorageOrder(StorageOrder order, StorageUpLo upLo,
81 IndexType n, IndexType k,
82 const MA *&A, IndexType &ldA)
83 {
84 IndexType kl = (upLo==Lower) ? k : 0;
85 IndexType ku = (upLo==Upper) ? k : 0;
86 switchBandStorageOrder(order, n, n, kl, ku, A, ldA);
87 }
88
89 template <typename IndexType, typename MA>
90 void
91 switchPackedStorageOrder(StorageOrder order, StorageUpLo &upLo,
92 IndexType n, const MA *&A)
93 {
94 assert(order==ColMajor);
95
96 MA *_A = new MA[n*(n+1)/2];
97 if (upLo==Upper) {
98 for (IndexType i=0; i<n; ++i) {
99 for (IndexType j=i; j<n; ++j) {
100 _A[j+i*(2*n-i-1)/2] = A[i+j*(j+1)/2];
101 }
102 }
103 } else {
104 for (IndexType i=0; i<n; ++i) {
105 for (IndexType j=0; j<=i; ++j) {
106 _A[j+i*(i+1)/2] = A[i+j*(2*n-j-1)/2];
107 }
108 }
109 }
110 A = _A;
111 }
112
113 template <typename IndexType, typename MA>
114 void
115 switchPackedStorageOrder(StorageOrder orderA, StorageUpLo &upLo,
116 IndexType n, const MA *A, MA *&_A)
117 {
118 if (upLo==Upper) {
119 if (orderA==ColMajor) {
120 for (IndexType i=0; i<n; ++i) {
121 for (IndexType j=i; j<n; ++j) {
122 _A[j+i*(2*n-i-1)/2] = A[i+j*(j+1)/2];
123 }
124 }
125 } else {
126 for (IndexType i=0; i<n; ++i) {
127 for (IndexType j=i; j<n; ++j) {
128 _A[i+j*(j+1)/2] = A[j+i*(2*n-i-1)/2];
129 }
130 }
131 }
132 } else {
133 if (orderA==ColMajor) {
134 for (IndexType i=0; i<n; ++i) {
135 for (IndexType j=0; j<=i; ++j) {
136 _A[j+i*(i+1)/2] = A[i+j*(2*n-j-1)/2];
137 }
138 }
139 } else {
140 for (IndexType i=0; i<n; ++i) {
141 for (IndexType j=0; j<=i; ++j) {
142 _A[i+j*(2*n-j-1)/2] = A[j+i*(i+1)/2];
143 }
144 }
145 }
146 }
147 }
148
149 template <typename IndexType, typename MA>
150 void
151 allocateFullStorage(StorageOrder order, IndexType m, IndexType n,
152 MA *&A, IndexType &ldA)
153 {
154 assert(order==RowMajor);
155
156 ldA = n;
157 A = new MA[ldA*m];
158 }
159
160 template <typename IndexType, typename MA>
161 void
162 allocatePackedStorage(IndexType n, MA *&A)
163 {
164 A = new MA[n*(n+1)/2];
165 }
166
167 template <typename MA>
168 void
169 releaseStorage(MA *A)
170 {
171 delete [] A;
172 }
173
174 #endif // CXXBLAS_INTERFACE_AUX_H