Lesson of 16 February 2010¶
Some simple C/C++ programs
Simple macro examples
File example6.cc
Example of FORTRAN linkage
File example7.cc
1/*
2 I am a comment
3 */
4
5// include the standard definition for I/O
6#include <iostream>
7// include the definition I/O manipulation command like setw, oct etc.
8#include <iomanip>
9
10// load the standard namespace
11using namespace std ;
12
13// define the fortran name as c name + underscore
14#define F77NAME(A) A##_
15
16typedef double valueType ;
17
18/* prototype of BLAS dgemm routine (fortran linkage) */
19
20extern "C" // telling the compiler that dgemm_ has C linkage convection
21void
22F77NAME(dgemm)( char const TRANSA[], // pass a string address
23 char const TRANSB[], // pass a string address
24 int const * M, // pass the address of the integer N
25 int const * N, // pass the address of the integer M
26 int const * K, // pass the address of the integer K
27 double const * ALPHA, // pass the address of the double ALPHA
28 double const A[], // pass the address of the vector/matrix A
29 int const * LDA, // pass the address of the row dimension of A
30 double const B[], // pass the address of the vector/matrix A
31 int const * LDB, // pass the address of the row dimension of B
32 double const * BETA, // pass the address of the double BETA
33 double C[], // pass the address of the row dimension of C
34 int const * LDC ) ; // pass the address of the vector/matrix A
35
36extern "C" // telling the compiler that sgemm_ has C linkage convection
37void
38F77NAME(sgemm)( char const TRANSA[],
39 char const TRANSB[],
40 int const * M,
41 int const * N,
42 int const * K,
43 float const * ALPHA,
44 float const A[],
45 int const * LDA,
46 float const B[],
47 int const * LDB,
48 float const * BETA,
49 float C[],
50 int const * LDC ) ;
51/*
52 FORTRAN do not support overloading.
53 We use inline + C++ overloading to simulate it.
54 In particular we define two routine "gemm" with
55 float and double argument respectively.
56 The use of inline results in a DIRECT call to the fortran
57 routine. No penalty in calling fotran routine!.
58 */
59
60inline
61void
62gemm( char const TRANSA[],
63 char const TRANSB[],
64 int const M,
65 int const N,
66 int const K,
67 float const ALPHA,
68 float const A[],
69 int const LDA,
70 float const B[],
71 int const LDB,
72 float const BETA,
73 float C[],
74 int const LDC ) {
75 F77NAME(sgemm)( TRANSA, TRANSB, &M, &N, &K, &ALPHA, A, &LDA, B, &LDB, &BETA, C, &LDC ) ;
76}
77
78inline
79void
80gemm( char const TRANSA[],
81 char const TRANSB[],
82 int const M,
83 int const N,
84 int const K,
85 double const ALPHA,
86 double const A[],
87 int const LDA,
88 double const B[],
89 int const LDB,
90 double const BETA,
91 double C[],
92 int const LDC ) {
93 F77NAME(dgemm)( TRANSA, TRANSB, &M, &N, &K, &ALPHA, A, &LDA, B, &LDB, &BETA, C, &LDC ) ;
94}
95
96int // the type of return value
97main() { // brackets contains the bosy of the program
98
99 // define matrix A and initialize
100 valueType A[3][3] = { {1, 2, 3},
101 {4, 5, 6},
102 {7, 8, 9} } ;
103
104 // define matrix B and initialize
105 valueType B[3][2] = { {1, 2},
106 {3, 4},
107 {5, 6} } ;
108
109 // define matrix C and DO NOT initialize
110 valueType C[10][10] ;
111
112 /* Due to the fact that C matrix are store by ROW and not by COLUMN
113 they are "seen" as transposed by the fortran routine.
114 Moreover (AB)^T = B^T A^T so that the
115 product of AB in C row major order is the
116 product BA in column major order!. */
117
118 // alpha * A * B + beta * C --> C
119 gemm( "N",
120 "N",
121 2,
122 3,
123 3,
124 1,
125 (valueType*)B, 2, // perform casting of the pointer to make compiler happy
126 (valueType*)A, 3, // perform casting of the pointer to make compiler happy
127 0,
128 (valueType*)C, // perform casting of the pointer to make compiler happy
129 10 ) ;
130
131 // display the result
132 for ( int i = 0 ; i < 3 ; ++i ) {
133 for ( int j = 0 ; j < 2 ; ++j ) {
134 cout << setw(10) << C[i][j] << " " ;
135 }
136 cout << '\n' ;
137 }
138
139 return 0 ; // return 0 to the OS
140}
141
142/*
143
144 Example of the assembler of the generated dgemm call
145
146 icc -S -fsource-asm -fcode-asm
147
148 ;;; gemm( "N",
149 ;;; "N",
150 ;;; 2,
151 ;;; 3,
152 ;;; 3,
153 ;;; 1,
154 ;;; (valueType*)B, 2, // perform casting of the pointer to make compiler happy
155 ;;; (valueType*)A, 3, // perform casting of the pointer to make compiler happy
156 ;;; 0,
157 ;;; (valueType*)C, // perform casting of the pointer to make compiler happy
158
159 movl $2, %ecx #example7.cc:119.3
160 movl $3, %r8d #example7.cc:119.3
161 lea 928(%rsp), %r9 #example7.cc:125.21
162 movq %r11, 952(%rsp) #example7.cc:105.21
163 lea 1088(%rsp), %rax #example7.cc:119.3
164 lea 976(%rsp), %r10 #example7.cc:126.21
165 lea 1096(%rsp), %r11 #example7.cc:119.3
166 movq %rdi, 960(%rsp) #example7.cc:105.21
167 lea 1056(%rsp), %rdi #example7.cc:119.3
168 movq %rsi, 968(%rsp) #example7.cc:105.21
169 lea 128(%rsp), %rsi #example7.cc:128.21
170 movl %ecx, 1064(%rsp) #example7.cc:119.3
171 movl %r8d, 1072(%rsp) #example7.cc:119.3
172 movl %r8d, 1080(%rsp) #example7.cc:119.3
173 movq %rdx, 1048(%rsp) #example7.cc:119.3
174 lea 1104(%rsp), %rdx #example7.cc:119.3
175 movl %ecx, 1088(%rsp) #example7.cc:119.3
176 lea 1072(%rsp), %rcx #example7.cc:119.3
177 movl %r8d, 1096(%rsp) #example7.cc:119.3
178 lea 1080(%rsp), %r8 #example7.cc:119.3
179 movq $0, 1056(%rsp) #example7.cc:119.3
180 movl $10, 1104(%rsp) #example7.cc:119.3
181 movq %r9, (%rsp) #example7.cc:125.21
182 lea 1048(%rsp), %r9 #example7.cc:119.3
183 movq %rax, 8(%rsp) #example7.cc:119.3
184 movq %r10, 16(%rsp) #example7.cc:126.21
185 movq %r11, 24(%rsp) #example7.cc:119.3
186 movq %rdi, 32(%rsp) #example7.cc:119.3
187 lea L_2__STRING.0(%rip), %rdi #example7.cc:119.3
188 movq %rsi, 40(%rsp) #example7.cc:128.21
189 movq %rdi, %rsi #example7.cc:119.3
190 movq %rdx, 48(%rsp) #example7.cc:119.3
191 lea 1064(%rsp), %rdx #example7.cc:119.3
192 call _dgemm_ #example7.cc:119.3
193
194*/
Some free numerical linear algebra library
The Lapack homepage
The Blas homepage
The ATLAS homepage