INGOR
Loading...
Searching...
No Matches
ytMatrix.h
1/*
2 math/ytMatrix.{h,c} : Matrix routines.
3 Copyright (C) 2018, Yoshinori Tamada <tamada A T ytlab.jp>
4 All rights reserved.
5
6 See LICENSE.txt for details of the licensing agreement.
7*/
8
9#ifndef __YTLIB_MATRIX_H
10#define __YTLIB_MATRIX_H
11
12/* Predefined constants */
13extern int ytMatrix_INT_ONE;
14extern int ytMatrix_INT_ZERO;
15extern double ytMatrix_DOUBLE_ONE;
16extern double ytMatrix_DOUBLE_ZERO;
17extern double ytMatrix_DOUBLE_MINUS_ONE;/* Minus One (-1)*/
18
19#include <stdio.h>
20
21void ytMatrix_dsyaxpy(const int n, double alpha, double * X, double * Y);
22double ytMatrix_trace(const int n, const int m, double * X);
23int ytMatrix_dsytrf_work_size(int n);
24int ytMatrix_dsyev_work_size(int n);
25void ytMatrix_print(FILE * fp, const double * X, int m, int n, int N);
26void ytMatrix_printTrans(FILE * fp, const double * X, int m, int n, int N);
27void ytMatrix_vprint(FILE * fp, const double * X, int n, int N);
28
29/* ------------------------------------------------------------------------- */
30/* Mac OS X vecLib */
31/* ------------------------------------------------------------------------- */
32#if defined(USE_MACOSX_VECLIB)
33
34#include <Accelerate/Accelerate.h>
35
36typedef __CLPK_integer ytMatrix_int;
37typedef __CLPK_doublereal ytMatrix_double;
38
39#ifndef USE_CBLAS_ROUTINES
40#define USE_CBLAS_ROUTINES
41#endif
42
43/* LAPACK Routines */
44extern char * ytMatrix_LU;
45extern char * ytMatrix_LL;
46extern char * ytMatrix_LN;
47extern char * ytMatrix_LV;
48extern char * ytMatrix_L_NO_TRANS;
49extern char * ytMatrix_L_TRANS;
50extern char * ytMatrix_LA;
51extern char * ytMatrix_LO;
52
53#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,LW,IF) dsytrf_(UL,N,A,LA,IP,W,LW,IF)
54#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) dsytri_(UL,N,A,LA,IP,W,IF)
55#define ytMatrix_dsytrs(UL,N,M,A,LA,IP,B,LB,IF) dsytrs_(UL,N,M,A,LA,IP,B,LB,IF)
56#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) dsyev_(JB,UL,N,A,LA,W,WO,LW,IF)
57#define ytMatrix_dgesv(N,M,A,LDA,IP,B,LDB,IF) dgesv_(N,M,A,LDA,IP,B,LDB,IF)
58#define ytMatrix_dgetrf(M,N,A,LA,IP,IF) dgetrf_(M,N,A,LA,IP,IF)
59#define ytMatrix_dgetri(N,A,LA,IP,W,LW,IF) dgetri_(N,A,LA,IP,W,LW,IF)
60#define ytMatrix_dgetrs(TA,N,M,A,LA,IP,B,LB,IF) dgetrs_(TA,N,M,A,LA,IP,B,LB,IF)
61#define ytMatrix_dgesvd(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF) \
62 dgesvd_(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF)
63
64/* USE_MACOSX_VECLIB */
65
66/* ------------------------------------------------------------------------- */
67/* INTEL MKL */
68/* ------------------------------------------------------------------------- */
69#elif defined(USE_INTEL_MKL)
70
71#include <mkl_blas.h>
72#include <mkl_lapack.h>
73#include <mkl_lapacke.h>
74
75typedef MKL_INT ytMatrix_int;
76typedef double ytMatrix_double;
77
78extern char * ytMatrix_T;
79extern char * ytMatrix_N;
80extern char * ytMatrix_C;
81extern char * ytMatrix_NO_TRANS;
82extern char * ytMatrix_TRANS;
83extern char * ytMatrix_U;
84extern char * ytMatrix_L;
85
86/* CBLAS ROUTINE */
87/* Level 1 BLAS */
88#define ytMatrix_dscal(N,ALPHA,X,IX) dscal(N,ALPHA,X,IX)
89#define ytMatrix_dcopy(N,X,IX,Y,IY) dcopy(N,X,IX,Y,IY)
90#define ytMatrix_daxpy(N,ALPHA,X,IX,Y,IY) daxpy(N,ALPHA,X,IX,Y,IY)
91#define ytMatrix_ddot(N,X,IX,Y,IY) ddot(N,X,IX,Y,IY)
92#define ytMatrix_dnrm2(N,X,IX) dnrm2(N,X,IX)
93#define ytMatrix_dasum(N,X,IX) dasum(N,X,IX)
94
95/* Level 2 BLAS */
96#define ytMatrix_dgemv(T,M,N,ALPHA,A,LDA,X,IX,BETA,Y,IY)\
97 dgemv(T,M,N,ALPHA,A,LDA,X,IX,BETA,Y,IY)
98#define ytMatrix_dsymv(UL,N,ALPHA,A,LA,X,IX,BETA,Y,IY)\
99 dsymv(UL,N,ALPHA,A,LA,X,IX,BETA,Y,IY)
100#define ytMatrix_dger(M,N,ALPHA,X,IX,Y,IY,A,LDA)\
101 dger(M,N,ALPHA,X,IX,Y,IY,A,LDA)
102
103/* Level 3 BLAS */
104#define ytMatrix_dgemm(TA,TB,M,N,K,ALPHA,A,LA,B,LB,BETA,C,LC) \
105 dgemm(TA,TB,M,N,K,ALPHA,A,LA,B,LB,BETA,C,LC)
106#define ytMatrix_dsyrk(UL,T,N,K,ALPHA,A,LDA,BETA,C,LC) \
107 dsyrk(UL,T,N,K,ALPHA,A,LDA,BETA,C,LC)
108
109/* LAPACK ROUTINE */
110extern char * ytMatrix_LU;
111extern char * ytMatrix_LL;
112extern char * ytMatrix_LN;
113extern char * ytMatrix_LV;
114extern char * ytMatrix_L_NO_TRANS;
115extern char * ytMatrix_L_TRANS;
116extern char * ytMatrix_LA;
117extern char * ytMatrix_LO;
118
119#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,LW,IF) \
120 LAPACKE_dsytrf(LAPACK_COL_MAJOR, *(UL), *(N), A, *(LA), IP)
121#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) \
122 LAPACKE_dsytri(LAPACK_COL_MAJOR, *(UL), *(N), A, *(LA), IP)
123//#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) dsytri_(UL,N,A,LA,IP,W,IF)
124//#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF)\
125// dsyev_(JB,UL,N,A,LA,W,WO,LW,IF)
126#define ytMatrix_dsytrs(UL,N,M,A,LA,IP,B,LB,IF) \
127 LAPACKE_dsytrs(LAPACK_COL_MAJOR,*(UL),*(N),*(M),A,*(LA),IP,B,*(LB))
128#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) \
129 LAPACKE_dsyev(LAPACK_COL_MAJOR,*(JB),*(UL),*(N),A,*(LA),W)
130#define ytMatrix_dgesv(N,M,A,LDA,IP,B,LDB,IF) dgesv_(N,M,A,LDA,IP,B,LDB,IF)
131#define ytMatrix_dgetrf(M,N,A,LA,IP,IF) dgetrf_(M,N,A,LA,IP,IF)
132#define ytMatrix_dgetri(N,A,LA,IP,W,LW,IF) dgetri_(N,A,LA,IP,W,LW,IF)
133#define ytMatrix_dgetrs(TA,N,M,A,LA,IP,B,LB,IF) dgetrs_(TA,N,M,A,LA,IP,B,LB,IF)
134#define ytMatrix_dgesvd(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF) \
135 dgesvd_(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF)
136
137
138/* ------------------------------------------------------------------------- */
139/* FUJITSU FCC/SSL */
140/* ------------------------------------------------------------------------- */
141#elif defined(USE_FUJITSU_FCC_SSL)
142
143#include <fj_lapack.h>
144
145typedef int ytMatrix_int;
146typedef double ytMatrix_double;
147
148extern char * ytMatrix_T;
149extern char * ytMatrix_N;
150extern char * ytMatrix_C;
151extern char * ytMatrix_NO_TRANS;
152extern char * ytMatrix_TRANS;
153extern char * ytMatrix_U;
154extern char * ytMatrix_L;
155extern char * ytMatrix_LEFT;
156extern char * ytMatrix_RIGHT;
157
158/* BLAS Routines */
159/* Level 1 BLAS */
160#define ytMatrix_dscal(N,ALPHA,X,IX) dscal_(N,ALPHA,X,IX)
161#define ytMatrix_dcopy(N,X,IX,Y,IY) dcopy_(N,X,IX,Y,IY)
162#define ytMatrix_daxpy(N,ALPHA,X,IX,Y,IY) daxpy_(N,ALPHA,X,IX,Y,IY)
163#define ytMatrix_ddot(N,X,IX,Y,IY) ddot_(N,X,IX,Y,IY)
164#define ytMatrix_dnrm2(N,X,IX) dnrm2_(N,X,IX)
165#define ytMatrix_dasum(N,X,IX) dasum_(N,X,IX)
166
167/* Level 2 BLAS */
168#define ytMatrix_dgemv(T,M,N,ALPHA,A,LDA,X,IX,BETA,Y,IY)\
169 dgemv_(T,M,N,ALPHA,A,LDA,X,IX,BETA,Y,IY,1)
170#define ytMatrix_dsymv(UL,N,ALPHA,A,LA,X,IX,BETA,Y,IY)\
171 dsymv_(UL,N,ALPHA,A,LA,X,IX,BETA,Y,IY,1)
172#define ytMatrix_dger(M,N,ALPHA,X,IX,Y,IY,A,LDA)\
173 dger_(M,N,ALPHA,X,IX,Y,IY,A,LDA)
174
175/* Level 3 BLAS */
176
177#define ytMatrix_dgemm(TA,TB,M,N,K,ALPHA,A,LA,B,LB,BETA,C,LC) \
178 dgemm_(TA,TB,M,N,K,ALPHA,A,LA,B,LB,BETA,C,LC,1,1)
179#define ytMatrix_dsymm(LR,UL,M,N,ALPHA,A,LA,B,LB,BETA,C,LC) \
180 dsymm_(LR,UL,M,N,ALPHA,A,LA,B,LB,BETA,C,LC,1,1)
181#define ytMatrix_dsyrk(UL,T,N,K,ALPHA,A,LDA,BETA,C,LC) \
182 dsyrk_(UL,T,N,K,ALPHA,A,LDA,BETA,C,LC,1,1)
183
184/* LAPACK Routines */
185extern char * ytMatrix_LU;
186extern char * ytMatrix_LL;
187extern char * ytMatrix_LN;
188extern char * ytMatrix_LV;
189extern char * ytMatrix_L_NO_TRANS;
190extern char * ytMatrix_L_TRANS;
191extern char * ytMatrix_LA;
192extern char * ytMatrix_LO;
193
194#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,IW,IF) dsytrf_(UL,N,A,LA,IP,W,IW,IF,1)
195#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) dsytri_(UL,N,A,LA,IP,W,IF,1)
196#define ytMatrix_dsytrs(UL,N,M,A,LA,IP,B,LB,IF) dsytrs_(UL,N,M,A,LA,IP,B,LB,IF,1)
197#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) dsyev_(JB,UL,N,A,LA,W,WO,LW,IF,1,1)
198#define ytMatrix_dgesv(N,M,A,LDA,IP,B,LDB,IF) dgesv_(N,M,A,LDA,IP,B,LDB,IF)
199#define ytMatrix_dgetrf(M,N,A,LA,IP,IF) dgetrf_(M,N,A,LA,IP,IF)
200#define ytMatrix_dgetri(N,A,LA,IP,W,LW,IF) dgetri_(N,A,LA,IP,W,LW,IF)
201#define ytMatrix_dgetrs(TA,N,M,A,LA,IP,B,LB,IF) dgetrs_(TA,N,M,A,LA,IP,B,LB,IF,1)
202#define ytMatrix_dgesvd(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF) \
203 dgesvd_(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF,1,1)
204
205/* End of FUJITSU FCC/SSL */
206
207
208#endif
209
210/* -- Arm Performace Libraries --------------------------------------------- */
211#if defined(USE_ARMPL)
212#include <armpl.h>
213#ifndef USE_CBLAS
214#define USE_CBLAS
215#endif
216
217extern char * ytMatrix_LU;
218extern char * ytMatrix_LL;
219extern char * ytMatrix_LN;
220extern char * ytMatrix_LV;
221extern char * ytMatrix_L_NO_TRANS;
222extern char * ytMatrix_L_TRANS;
223extern char * ytMatrix_LA;
224extern char * ytMatrix_LO;
225
226#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,LW,IF) dsytrf_(UL,N,A,LA,IP,W,LW,IF)
227#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) dsytri_(UL,N,A,LA,IP,W,IF)
228#define ytMatrix_dsytrs(UL,N,M,A,LA,IP,B,LB,IF) dsytrs_(UL,N,M,A,LA,IP,B,LB,IF)
229#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) dsyev_(JB,UL,N,A,LA,W,WO,LW,IF)
230#define ytMatrix_dgesv(N,M,A,LDA,IP,B,LDB,IF) dgesv_(N,M,A,LDA,IP,B,LDB,IF)
231#define ytMatrix_dgetrf(M,N,A,LA,IP,IF) dgetrf_(M,N,A,LA,IP,IF)
232#define ytMatrix_dgetri(N,A,LA,IP,W,LW,IF) dgetri_(N,A,LA,IP,W,LW,IF)
233#define ytMatrix_dgetrs(TA,N,M,A,LA,IP,B,LB,IF) dgetrs_(TA,N,M,A,LA,IP,B,LB,IF)
234#define ytMatrix_dgesvd(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF) \
235 dgesvd_(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF)
236
237#endif
238
239/* -- CBLAS ---------------------------------------------------------------- */
240#if defined(USE_CBLAS)
241#include <cblas.h>
242#ifndef USE_CBLAS_ROUTINES
243#define USE_CBLAS_ROUTINES
244#endif /* ifndef USE_CBLAS_ROUTINES */
245#endif /* USE_CBLAS */
246
247#if defined(USE_CBLAS_ROUTINES)
248
249#define ytMatrix_T CblasTrans
250#define ytMatrix_N CblasNoTrans
251#define ytMatrix_C CblasConjTrans
252
253#define ytMatrix_NO_TRANS CblasNoTrans
254#define ytMatrix_TRANS CblasTrans
255
256#define ytMatrix_U CblasUpper
257#define ytMatrix_L CblasLower
258
259#define ytMatrix_LEFT CblasLeft
260#define ytMatrix_RIGHT CblasRight
261
262#define ytMatrix_dscal(N,ALPHA,X,IX) cblas_dscal(*N,*ALPHA,X,*IX)
263#define ytMatrix_dcopy(N,X,IX,Y,IY) cblas_dcopy(*N,X,*IX,Y,*IY)
264#define ytMatrix_daxpy(N,ALPHA,X,IX,Y,IY) cblas_daxpy(*N,*ALPHA,X,*IX,Y,*IY)
265#define ytMatrix_ddot(N,X,IX,Y,IY) cblas_ddot(*N,X,*IX,Y,*IY)
266#define ytMatrix_dnrm2(N,X,IX) cblas_dnrm2(*N,X,*IX)
267#define ytMatrix_dasum(N,X,IX) cblas_dasum(*N,X,*IX)
268
269/* Level 2 BLAS */
270#define ytMatrix_dgemv(T,M,N,ALPHA,A,LDA,X,IX,BETA,Y,IY)\
271 cblas_dgemv(CblasColMajor,T,*M,*N,*ALPHA,A,*LDA,X,*IX,*BETA,Y,*IY)
272#define ytMatrix_dsymv(UL,N,ALPHA,A,LA,X,IX,BETA,Y,IY)\
273 cblas_dsymv(CblasColMajor,UL,*N,*ALPHA,A,*LA,X,*IX,*BETA,Y,*IY)
274#define ytMatrix_dger(M,N,ALPHA,X,IX,Y,IY,A,LDA)\
275 cblas_dger(CblasColMajor,*M,*N,*ALPHA,X,*IX,Y,*IY,A,*LDA)
276
277/* Level 3 BLAS */
278#define ytMatrix_dgemm(TA,TB,M,N,K,ALPHA,A,LA,B,LB,BETA,C,LC) \
279 cblas_dgemm(CblasColMajor,TA,TB,*M,*N,*K,*ALPHA,A,*LA,B,*LB,*BETA,C,*LC)
280#define ytMatrix_dsymm(LR,UL,M,N,ALPHA,A,LA,B,LB,BETA,C,LC) \
281 cblas_dsymm(CblasColMajor,LR,UL,*M,*N,*ALPHA,A,*LA,B,*LB,*BETA,C,*LC)
282#define ytMatrix_dsyrk(UL,T,N,K,ALPHA,A,LDA,BETA,C,LC) \
283 cblas_dsyrk(CblasColMajor,UL,T,*N,*K,*ALPHA,A,*LDA,*BETA,C,*LC)
284
285#endif /* USE_CBLAS_ROUTINES */
286
287#ifdef USE_LAPACKE
288
289#include <lapacke.h>
290extern char * ytMatrix_LU;
291extern char * ytMatrix_LL;
292extern char * ytMatrix_LN;
293extern char * ytMatrix_LV;
294extern char * ytMatrix_L_NO_TRANS;
295extern char * ytMatrix_L_TRANS;
296extern char * ytMatrix_LA;
297extern char * ytMatrix_LO;
298
299/* LAPACK ROUTINES */
300#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,LW,IF) \
301 LAPACKE_dsytrf(LAPACK_COL_MAJOR, *(UL), *(N), A, *(LA), IP)
302#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) \
303 LAPACKE_dsytri(LAPACK_COL_MAJOR, *(UL), *(N), A, *(LA), IP)
304#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) \
305 LAPACKE_dsyev(LAPACK_COL_MAJOR,*(JB),*(UL),*(N),A,*(LA),W)
306#endif /* End of USE_LAPACKE */
307
308/* ------------------------------------------------------------------------- */
309
310#ifdef USE_LAPACK
311/* LAPACK Routines */
312extern char * ytMatrix_LU;
313extern char * ytMatrix_LL;
314extern char * ytMatrix_LN;
315extern char * ytMatrix_LV;
316extern char * ytMatrix_L_NO_TRANS;
317extern char * ytMatrix_L_TRANS;
318extern char * ytMatrix_LA;
319extern char * ytMatrix_LO;
320
321void dsytrf_(char * UL, int * N, double * A, int * ldA, int * IP, double * WORK, int * LWORK, int * info);
322void dsytri_(char * UL, int * N, double * A, int * ldA, int * IP, double * W, int * info);
323void dsyev_(char * JB, char * UL, int * N, double * A, int * ldA, double * W,
324 double * WORK, int * LWORK, int * info);
325void dgesv_(int * N, int * M, double * A, int * ldA, int * IP, double * B, int * ldB, int * info);
326void dgetrf_(int * M, int * N, double * A, int * ldA, int * IP, int * info);
327void dgetri_(int * N, double * A, int * ldA, int * IP, double * WORK, int * LWORK, int * info);
328void dgetrs_(char * TA, int * N, int * M, double * A, int * ldA, int * IP,
329 double * B, int * ldB, int * info);
330void dgesvd_(char * JU, char * JV, int * M, int * N, double * A, int * ldA, double * S,
331 double * U, int * ldU, double * VT, int * ldVT, double * WORK, int * LWORK, int * info);
332
333#define ytMatrix_dsytrf(UL,N,A,LA,IP,W,LW,IF) dsytrf_(UL,N,A,LA,IP,W,LW,IF)
334#define ytMatrix_dsytri(UL,N,A,LA,IP,W,IF) dsytri_(UL,N,A,LA,IP,W,IF)
335#define ytMatrix_dsytrs(UL,N,M,A,LA,IP,B,LB,IF) dsytrs_(UL,N,M,A,LA,IP,B,LB,IF)
336#define ytMatrix_dsyev(JB,UL,N,A,LA,W,WO,LW,IF) dsyev_(JB,UL,N,A,LA,W,WO,LW,IF)
337#define ytMatrix_dgesv(N,M,A,LDA,IP,B,LDB,IF) dgesv_(N,M,A,LDA,IP,B,LDB,IF)
338#define ytMatrix_dgetrf(M,N,A,LA,IP,IF) dgetrf_(M,N,A,LA,IP,IF)
339#define ytMatrix_dgetri(N,A,LA,IP,W,LW,IF) dgetri_(N,A,LA,IP,W,LW,IF)
340#define ytMatrix_dgetrs(TA,N,M,A,LA,IP,B,LB,IF) dgetrs_(TA,N,M,A,LA,IP,B,LB,IF)
341#define ytMatrix_dgesvd(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF) \
342 dgesvd_(JU,JV,M,N,A,LA,S,U,LU,VT,LVT,W,LW,IF)
343
344#endif /* USE_LAPACK */
345
346
347
348#endif /* __YTLIB_MATRIX_H */
ytMatrix_dsytrf_work_size
int ytMatrix_dsytrf_work_size(int n)
Returns the optimal working memory size for dsytrf routine.
Definition ytMatrix.c:106
ytMatrix_trace
double ytMatrix_trace(const int n, const int m, double *X)
Returns the trace of the given matrix.
Definition ytMatrix.c:89
ytMatrix_print
void ytMatrix_print(FILE *fp, const double *X, int m, int n, int N)
Prints the matrix.
Definition ytMatrix.c:151
ytMatrix_printTrans
void ytMatrix_printTrans(FILE *fp, const double *X, int m, int n, int N)
Prints the transposed or row-major matrix.
Definition ytMatrix.c:168
ytMatrix_L_NO_TRANS
@ ytMatrix_L_NO_TRANS
No transpose.
Definition ytMatrix.c:63
ytMatrix_L_TRANS
@ ytMatrix_L_TRANS
Transpose.
Definition ytMatrix.c:66
Generated by doxygen 1.9.8