blasLapackWrapper.h 6.36 KB
Newer Older
Baum's avatar
Baum committed
1
// Xerus - A General Purpose Tensor Library
Sebastian Wolf's avatar
Sebastian Wolf committed
2
// Copyright (C) 2014-2016 Benjamin Huber and Sebastian Wolf. 
Baum's avatar
Baum committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

20
/**
21
22
23
* @file
* @brief Header file for the blas and lapack wrapper functions.
*/
24

Baum's avatar
Baum committed
25
26
#pragma once

27
#include "misc/standard.h"
28
#include <memory>
Baum's avatar
Baum committed
29
30


31
namespace xerus {
Sebastian Wolf's avatar
Sebastian Wolf committed
32
	/**
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
	* @brief In this namespace the minimal wrappers for the BLAS and LAPACK functions are collected.
	* @details As an end user of xerus it should never be nessecary to call any of these functions, unless
	* a seriously low level implementation of a critical part of an algorithm is required.
	*/
	namespace blasWrapper {
		
		//----------------------------------------------- LEVEL I BLAS ----------------------------------------------------------
		
		///@brief: Computes the two norm =||x||
		double two_norm(const double* const _x, const size_t _n);
		
		///@brief: Computes the dot product = x^T*y
		double dot_product(const double* const _x, const size_t _n, const double* const _y);
		
		
		//----------------------------------------------- LEVEL II BLAS ---------------------------------------------------------
		
		///@brief: Perfroms x = alpha*OP(A)*y
		void matrix_vector_product(double* const _x, const size_t _m, const double _alpha, const double* const _A, const size_t _n, const bool _transposed, const double* const _y);
		
		///@brief: Performs A = alpha*x*y^T
		void dyadic_vector_product(double* _A, const size_t _m, const size_t _n, const double _alpha, const double*const _x, const double*const _y);
		
		//----------------------------------------------- LEVEL III BLAS --------------------------------------------------------
		///@brief: Performs the Matrix-Matrix product C = alpha*OP(A) * OP(B)
		void matrix_matrix_product( double* const _C,
									const size_t _leftDim,
									const size_t _rightDim,
									const double _alpha,
									const double* const _A,
									const size_t _lda,
									const bool _transposeA,
									const size_t _middleDim,
									const double* const _B,
									const size_t _ldb,
									const bool _transposeB);
		
		///@brief: Performs the Matrix-Matrix product C = alpha*OP(A) * OP(B)
71
		static XERUS_force_inline void matrix_matrix_product( double* const _C,
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
									const size_t _leftDim,
									const size_t _rightDim,
									const double _alpha,
									const double* const _A,
									const bool _transposeA,
									const size_t _middleDim,
									const double* const _B,
									const bool _transposeB) {
			
			matrix_matrix_product( _C, _leftDim, _rightDim, _alpha, _A, _transposeA ? _leftDim : _middleDim, _transposeA, _middleDim, _B, _transposeB ? _middleDim : _rightDim, _transposeB);
		}
		
		//----------------------------------------------- LAPACK ----------------------------------------------------------------
		
		///@brief: Performs (U,S,V) = SVD(A)
		void svd( double* const _U, double* const _S, double* const _Vt, const double* const _A, const size_t _m, const size_t _n);
Baum's avatar
Baum committed
88

89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
		///@brief: Performs (U,S,V) = SVD(A). Destroys A.
		void svd_destructive( double* const _U, double* const _S, double* const _Vt, double* const _A, const size_t _m, const size_t _n);
		
		
		///@brief: splits A = Q*C, with @a _C an rxn matrix (where r is the rank of @a _A) and @a _Q orthogonal.
		std::tuple<std::unique_ptr<double[]>, std::unique_ptr<double[]>, size_t> qc(const double* const _A, const size_t _m, const size_t _n);
		
		///@brief: splits A = Q*C, with @a _C an rxn matrix (where r is the rank of @a _A) and @a _Q orthogonal. Destroys A.
		std::tuple<std::unique_ptr<double[]>, std::unique_ptr<double[]>, size_t> qc_destructive(double* const _A, const size_t _m, const size_t _n);
		
		
		///@brief: splits A = C*Q, with @a _C an rxm matrix (where r is the rank of @a _A) and @a _Q orthogonal.
		std::tuple<std::unique_ptr<double[]>, std::unique_ptr<double[]>, size_t> cq(const double* const _A, const size_t _m, const size_t _n);
		
		///@brief: splits A = C*Q, with @a _C an rxm matrix (where r is the rank of @a _A) and @a _Q orthogonal. Destroys A.
		std::tuple<std::unique_ptr<double[]>, std::unique_ptr<double[]>, size_t> cq_destructive(double* const _A, const size_t _m, const size_t _n);
		
		
		///@brief: Performs (Q,R) = QR(A)
		void qr( double* const _Q, double* const _R, const double* const _A, const size_t _m, const size_t _n);
		
		///@brief: Performs (AtoQ,R) = QR(AtoQ)
		void inplace_qr(double* const _AtoQ, double* const _R, const size_t _m, const size_t _n);
		
		///@brief: Performs (Q,R) = QR(A), destroys A in the process
		void qr_destructive( double* const _Q, double* const _R, double* const _A, const size_t _m, const size_t _n);
		
		
		
		///@brief: Performs (R,Q) = RQ(A)
		void rq( double* const _R, double* const _Q, const double* const _A, const size_t _m, const size_t _n);
		
		///@brief: Performs (R,AtoQ) = RQ(AtoQ)
		void inplace_rq( double* const _R, double* const _AtoQ, const size_t _m, const size_t _n);
		
		///@brief: Performs (R,Q) = RQ(A), destroys A in the process
		void rq_destructive( double* const _R, double* const _Q, double* const _A, const size_t _m, const size_t _n);
126

127
128
		
		///@brief: Solves Ax = b for x
129
		void solve(double* _x, const double* _A, size_t _m, size_t _n, const double* _b, size_t _p);
130
131
132
133
		
		
		
		///@brief: Solves min ||Ax - b||_2 for x
134
		void solve_least_squares( double* const _x, const double* const _A, const size_t _m, const size_t _n, const double* const _b, const size_t _p);
135
136
		
		///@brief: Solves min ||Ax - b||_2 for x
137
		void solve_least_squares_destructive( double* const _x, double* const _A, const size_t _m, const size_t _n, double* const _b, const size_t _p);
138
139
140
141
142
143
		
		
		
		///@brief: Performs B = A^-1
		void inverse( double* const _B, const double* const _A, const size_t _m, const size_t _n);
	}
Baum's avatar
Baum committed
144

145
}