Lesson of 16 July 2008¶
Solving a Boundary Value Problem with the use of Lapack routines
Solve a simple BVP problem
File bvp/main.cc
:open:
1/*
2 A simple program that compute approximate solution of
3 the boundary value problem
4
5 y'' + p y' + q y = r
6 y(a) = ya
7 y(b) = yb
8
9 by finite difference
10
11*/
12
13#include <iostream>
14#include <fstream>
15#include <cmath>
16
17// from STL
18#include <string>
19#include <vector>
20
21using namespace std ;
22
23typedef double valueType ;
24typedef int indexType ;
25
26/*
27 Parameter for the specific problem
28
29 y'' + y' - x*y = -x^2 ;
30
31 */
32
33static valueType p( valueType x ) { return 1 ; }
34static valueType q( valueType x ) { return -x ; }
35static valueType r( valueType x ) { return -x*x ; }
36
37/*
38 Build the coefficient of the tridiagonal matrix
39 for the problem
40
41 // first row
42 alpha[0] * y[0] + gamma[0] * y[1] = omega[0] - beta[0] * ya
43
44 // internal rows
45 beta[k-1] * y[k-1] + alpha[k] * y[k] + gamma[k] * y[k+1] = omega[k]
46
47 // last row
48 beta[n-2] * y[n-2] + alpha[n-1] * y[n-1] = omega[n-1] - gamma[n-1] * yb
49
50 alpha[k] =
51
52
53 | alpha[0] gamma[0] 0 0 0 |
54 | beta[0] alpha[1] gamma[1] 0 0 |
55 T = | 0 beta[1] alpha[2] gamma[2] 0 |
56 | 0 0 beta[2] alpha[3] gamma[3] |
57 | 0 0 0 beta[3] alpha[4] |
58
59
60 b = [ omega[0], omega[1],...., omega[n-2] ]^T ;
61
62 alpha[k] = -2/h^2 + q(x[k+1])
63 beta[k] = 1/h^2 + p(x[k+2])/(2h)
64 gamma[k] = 1/h^2 - p(x[k+1])/(2h)
65 omega[k] = r(x[k+1])
66 omega[0] = r(x[1])-beta[0]*ya
67 omega[N-2] = r(x[N-1])-gamma[N-2]*yb
68
69 */
70
71static
72void
73buildSystem( indexType const N, // number of interval of subdivision
74 valueType const a, // left boundary
75 valueType const ya, // value of the solution on left boundary
76 valueType const b, // right boundary
77 valueType const yb, // value of the solution on right boundary
78 vector<valueType> & alpha, // main diagonal
79 vector<valueType> & beta, // lower diagonal
80 vector<valueType> & gamma, // upper diagonal
81 vector<valueType> & omega // r.h.s.
82 ) {
83 // some check
84 if ( b <= a ) throw "bad interval" ; // interrupt the code with an exception
85 if ( N <= 1 ) throw "N must be > 0" ;
86
87 // compute the cell size
88 valueType h = (b-a)/N ;
89
90 // initialize the vectors (to size 0)
91 alpha . clear() ;
92 beta . clear() ;
93 gamma . clear() ;
94 omega . clear() ;
95
96 // reserve memory to avoid multiple allocation
97 alpha . reserve(N-1) ;
98 beta . reserve(N-2) ;
99 gamma . reserve(N-2) ;
100 omega . reserve(N-1) ;
101
102 // fill the vector alpha and omega
103 for ( indexType k = 0 ; k < N-1 ; ++k ) {
104 valueType xkp1 = a + h * (k+1) ;
105 alpha . push_back( -2/(h*h) + q(xkp1) ) ;
106 omega . push_back( r(xkp1) ) ;
107 }
108 // fill the vector beta
109 for ( indexType k = 0 ; k < N-2 ; ++k ) {
110 valueType xkp2 = a + h * (k+2) ;
111 beta . push_back( 1/(h*h) + p(xkp2)/(2*h) ) ;
112 }
113 // fill the vector gamma
114 for ( indexType k = 0 ; k < N-2 ; ++k ) {
115 valueType xkp1 = a + h * (k+1) ;
116 gamma . push_back( 1/(h*h) - p(xkp1)/(2*h) ) ;
117 }
118
119 // setup boundary conditions
120 omega.front() -= beta.front() * ya ;
121 omega.back() -= gamma.back() * yb ;
122}
123
124/*
125
126 solve the tridiagonal system
127
128
129 | alpha[0] gamma[0] 0 0 0 |
130 | beta[0] alpha[1] gamma[1] 0 0 |
131 T = | 0 beta[1] alpha[2] gamma[2] 0 |
132 | 0 0 beta[2] alpha[3] gamma[3] |
133 | 0 0 0 beta[3] alpha[4] |
134
135
136 b = [ omega[0], omega[1],...., omega[n-2] ]^T ;
137
138 */
139
140/*
141 Access to the FORTRAN routine DGTSV
142
143 */
144
145extern "C" void dgtsv( int * N,
146 int * NRHS,
147 double * DL,
148 double * D,
149 double * DU,
150 double * B,
151 int * LDB,
152 int * INFO );
153
154static
155void
156solveSystem( vector<valueType> const & alpha, // main diagonal
157 vector<valueType> const & beta, // lower diagonal
158 vector<valueType> const & gamma, // upper diagonal
159 vector<valueType> const & omega, // r.h.s.
160 vector<valueType> & sol // the solution
161 ) {
162 // some check
163 indexType N = alpha.size()+1 ;
164 if ( beta.size() != N-2 ) throw "bad size for beta " ; // interrupt the code with an exception
165 if ( gamma.size() != N-2 ) throw "bad size for gamma " ; // interrupt the code with an exception
166 if ( omega.size() != N-1 ) throw "bad size for omega " ; // interrupt the code with an exception
167
168 sol . resize( N-1 ) ;
169 copy( omega.begin(), omega.end(), sol.begin() ) ;
170
171 // working vector
172 vector<valueType> D(alpha), DU(gamma), DL(beta) ;
173
174 indexType n = N-1 ;
175 indexType one = 1 ;
176 indexType INFO ;
177 dgtsv( &n, &one, &DL.front(),
178 &D.front(),
179 &DU.front(),
180 &sol.front(),
181 &n,
182 &INFO ) ;
183
184 /* INFO (output) INTEGER
185 * = 0: successful exit
186 * < 0: if INFO = -i, the i-th argument had an illegal value
187 * > 0: if INFO = i, U(i,i) is exactly zero, and the solution
188 * has not been computed. The factorization has not been
189 * completed unless i = N.
190 */
191 if ( INFO < 0 ) throw "bad argument for dgtsv" ;
192 if ( INFO > 0 ) throw "singular tridiagonal matrix for for dgtsv" ;
193
194}
195
196
197int
198main( int argc, char * argv[] ) {
199 indexType N ;
200
201 if ( argc > 2 ) {
202 // too nuch argument, issue an error
203 cerr << "usage: " << argv[0] << " N\n" ;
204 return 1 ;
205 }
206
207 if ( argc == 1 ) {
208 // only one argument, ask to the user the file name
209 cout << "Number of interval: " ;
210 cin >> N ;
211 } else {
212 // two argument, get file name from command line
213 scanf( argv[1], "%d", &N ) ;
214 }
215
216 // open the file
217 ofstream file( "result.txt" ) ;
218
219 // check if the file exist is opened etc.
220 if ( file . bad() ) {
221 cerr << "error in openeing : result.txt\n" ;
222 return 1 ;
223 }
224
225 valueType a = 0 ;
226 valueType ya = 0 ;
227 valueType b = 1 ;
228 valueType yb = 2 ;
229 vector<valueType> alpha, beta, gamma, omega, sol ;
230
231 try {
232 buildSystem( N, a, ya, b, yb, alpha, beta, gamma, omega ) ;
233 solveSystem( alpha, beta, gamma, omega, sol ) ;
234 // add initial and final boundary value
235 sol . push_back( yb ) ;
236 sol . insert( sol.begin(), ya ) ;
237 }
238 catch ( char const msg[] ) {
239 cerr << "Error found: " << msg << "\n" ;
240 exit(0) ;
241 }
242 catch ( ... ) { // catch any kind of error
243 cerr << "Unknown error\n" ;
244 exit(0) ;
245 }
246 file << "x\ty\n" ;
247 indexType k = 0 ;
248 for ( vector<valueType>::const_iterator i = sol.begin() ;
249 i != sol.end() ; ++i )
250 file << (a + (k++)*(b-a)/N) << "\t" << *i << "\n" ;
251
252 // close the file
253 file . close() ;
254}
Solve a simple BVP problem with a more customizable interface
File bvpbetter/main.cc
:open:
1/*
2 A simple program that compute approximate solution of
3 the boundary value problem
4
5 y'' + p y' + q y = r
6 y(a) = ya
7 y(b) = yb
8
9 by finite difference
10
11*/
12
13#include <iostream>
14#include <fstream>
15#include <cmath>
16
17// from STL
18#include <string>
19#include <vector>
20
21#include "Blas.hh"
22#include "Lapack.hh"
23
24using namespace std ;
25
26typedef double valueType ;
27typedef int indexType ;
28
29/*
30 Parameter for the specific problem
31
32 y'' + y' - x*y = -x^2 ;
33
34 */
35
36static valueType p( valueType x ) { return 1 ; }
37static valueType q( valueType x ) { return -x ; }
38static valueType r( valueType x ) { return -x*x ; }
39
40/*
41 Build the coefficient of the tridiagonal matrix
42 for the problem
43
44 // first row
45 alpha[0] * y[0] + gamma[0] * y[1] = omega[0] - beta[0] * ya
46
47 // internal rows
48 beta[k-1] * y[k-1] + alpha[k] * y[k] + gamma[k] * y[k+1] = omega[k]
49
50 // last row
51 beta[n-2] * y[n-2] + alpha[n-1] * y[n-1] = omega[n-1] - gamma[n-1] * yb
52
53 alpha[k] =
54
55
56 | alpha[0] gamma[0] 0 0 0 |
57 | beta[0] alpha[1] gamma[1] 0 0 |
58 T = | 0 beta[1] alpha[2] gamma[2] 0 |
59 | 0 0 beta[2] alpha[3] gamma[3] |
60 | 0 0 0 beta[3] alpha[4] |
61
62
63 b = [ omega[0], omega[1],...., omega[n-2] ]^T ;
64
65 alpha[k] = -2/h^2 + q(x[k+1])
66 beta[k] = 1/h^2 + p(x[k+2])/(2h)
67 gamma[k] = 1/h^2 - p(x[k+1])/(2h)
68 omega[k] = r(x[k+1])
69 omega[0] = r(x[1])-beta[0]*ya
70 omega[N-2] = r(x[N-1])-gamma[N-2]*yb
71
72 */
73
74static
75void
76buildSystem( indexType const N, // number of interval of subdivision
77 valueType const a, // left boundary
78 valueType const ya, // value of the solution on left boundary
79 valueType const b, // right boundary
80 valueType const yb, // value of the solution on right boundary
81 vector<valueType> & alpha, // main diagonal
82 vector<valueType> & beta, // lower diagonal
83 vector<valueType> & gamma, // upper diagonal
84 vector<valueType> & omega // r.h.s.
85 ) {
86 // some check
87 if ( b <= a ) throw "bad interval" ; // interrupt the code with an exception
88 if ( N <= 1 ) throw "N must be > 0" ;
89
90 // compute the cell size
91 valueType h = (b-a)/N ;
92
93 // initialize the vectors (to size 0)
94 alpha . clear() ;
95 beta . clear() ;
96 gamma . clear() ;
97 omega . clear() ;
98
99 // reserve memory to avoid multiple allocation
100 alpha . reserve(N-1) ;
101 beta . reserve(N-2) ;
102 gamma . reserve(N-2) ;
103 omega . reserve(N-1) ;
104
105 // fill the vector alpha and omega
106 for ( indexType k = 0 ; k < N-1 ; ++k ) {
107 valueType xkp1 = a + h * (k+1) ;
108 alpha . push_back( -2/(h*h) + q(xkp1) ) ;
109 omega . push_back( r(xkp1) ) ;
110 }
111 // fill the vector beta
112 for ( indexType k = 0 ; k < N-2 ; ++k ) {
113 valueType xkp2 = a + h * (k+2) ;
114 beta . push_back( 1/(h*h) + p(xkp2)/(2*h) ) ;
115 }
116 // fill the vector gamma
117 for ( indexType k = 0 ; k < N-2 ; ++k ) {
118 valueType xkp1 = a + h * (k+1) ;
119 gamma . push_back( 1/(h*h) - p(xkp1)/(2*h) ) ;
120 }
121
122 // setup boundary conditions
123 omega.front() -= beta.front() * ya ;
124 omega.back() -= gamma.back() * yb ;
125}
126
127/*
128
129 solve the tridiagonal system
130
131
132 | alpha[0] gamma[0] 0 0 0 |
133 | beta[0] alpha[1] gamma[1] 0 0 |
134 T = | 0 beta[1] alpha[2] gamma[2] 0 |
135 | 0 0 beta[2] alpha[3] gamma[3] |
136 | 0 0 0 beta[3] alpha[4] |
137
138
139 b = [ omega[0], omega[1],...., omega[n-2] ]^T ;
140
141 */
142
143/*
144 Access to the FORTRAN routine DGTSV
145
146 */
147
148extern "C" void sgtsv( int * N,
149 int * NRHS,
150 float * DL,
151 double * D,
152 double * DU,
153 double * B,
154 int * LDB,
155 int * INFO );
156
157static
158void
159solveSystem( vector<valueType> const & alpha, // main diagonal
160 vector<valueType> const & beta, // lower diagonal
161 vector<valueType> const & gamma, // upper diagonal
162 vector<valueType> const & omega, // r.h.s.
163 vector<valueType> & sol // the solution
164 ) {
165 // some check
166 indexType N = alpha.size()+1 ;
167 if ( beta.size() != N-2 ) throw "bad size for beta " ; // interrupt the code with an exception
168 if ( gamma.size() != N-2 ) throw "bad size for gamma " ; // interrupt the code with an exception
169 if ( omega.size() != N-1 ) throw "bad size for omega " ; // interrupt the code with an exception
170
171 sol . resize( N-1 ) ;
172 copy( omega.begin(), omega.end(), sol.begin() ) ;
173
174 // working vector
175 vector<valueType> D(alpha), DU(gamma), DL(beta) ;
176
177 indexType INFO = lapack::gtsv( N-1, 1,
178 &DL.front(),
179 &D.front(),
180 &DU.front(),
181 &sol.front(),
182 N-1 ) ;
183
184 /* INFO (output) INTEGER
185 * = 0: successful exit
186 * < 0: if INFO = -i, the i-th argument had an illegal value
187 * > 0: if INFO = i, U(i,i) is exactly zero, and the solution
188 * has not been computed. The factorization has not been
189 * completed unless i = N.
190 */
191 if ( INFO < 0 ) throw "bad argument for dgtsv" ;
192 if ( INFO > 0 ) throw "singular tridiagonal matrix for for dgtsv" ;
193
194}
195
196
197int
198main( int argc, char * argv[] ) {
199 indexType N ;
200
201 if ( argc > 2 ) {
202 // too nuch argument, issue an error
203 cerr << "usage: " << argv[0] << " N\n" ;
204 return 1 ;
205 }
206
207 if ( argc == 1 ) {
208 // only one argument, ask to the user the file name
209 cout << "Number of interval: " ;
210 cin >> N ;
211 } else {
212 // two argument, get file name from command line
213 scanf( argv[1], "%d", &N ) ;
214 }
215
216 // open the file
217 ofstream file( "result.txt" ) ;
218
219 // check if the file exist is opened etc.
220 if ( file . bad() ) {
221 cerr << "error in openeing : result.txt\n" ;
222 return 1 ;
223 }
224
225 valueType a = 0 ;
226 valueType ya = 0 ;
227 valueType b = 1 ;
228 valueType yb = 2 ;
229 vector<valueType> alpha, beta, gamma, omega, sol ;
230
231 try {
232 buildSystem( N, a, ya, b, yb, alpha, beta, gamma, omega ) ;
233 solveSystem( alpha, beta, gamma, omega, sol ) ;
234 // add initial and final boundary value
235 sol . push_back( yb ) ;
236 sol . insert( sol.begin(), ya ) ;
237 }
238 catch ( char const msg[] ) {
239 cerr << "Error found: " << msg << "\n" ;
240 exit(0) ;
241 }
242 catch ( ... ) { // catch any kind of error
243 cerr << "Unknown error\n" ;
244 exit(0) ;
245 }
246 file << "x\ty\n" ;
247 indexType k = 0 ;
248 for ( vector<valueType>::const_iterator i = sol.begin() ;
249 i != sol.end() ; ++i )
250 file << (a + (k++)*(b-a)/N) << "\t" << *i << "\n" ;
251
252 // close the file
253 file . close() ;
254}
Blas C/FORTRAN interaface
File bvpbetter/Blas.hh
:open:
1/*--------------------------------------------------------------------------*\
2 | |
3 | This program is free software; you can redistribute it and/or modify |
4 | it under the terms of the GNU General Public License as published by |
5 | the Free Software Foundation; either version 2, or (at your option) |
6 | any later version. |
7 | |
8 | This program is distributed in the hope that it will be useful, |
9 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
10 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
11 | GNU General Public License for more details. |
12 | |
13 | You should have received a copy of the GNU General Public License |
14 | along with this program; if not, write to the Free Software |
15 | Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
16 | |
17 | Copyright (C) 2008 |
18 | |
19 | , __ , __ |
20 | /|/ \ /|/ \ |
21 | | __/ _ ,_ | __/ _ ,_ |
22 | | \|/ / | | | | \|/ / | | | |
23 | |(__/|__/ |_/ \_/|/|(__/|__/ |_/ \_/|/ |
24 | /| /| |
25 | \| \| |
26 | |
27 | Enrico Bertolazzi |
28 | Dipartimento di Ingegneria Meccanica e Strutturale |
29 | Universita` degli Studi di Trento |
30 | Via Mesiano 77, I-38050 Trento, Italy |
31 | email: enrico.bertolazzi@unitn.it |
32 | |
33\*--------------------------------------------------------------------------*/
34
35///
36/// file: Blas.hh
37///
38
39#ifndef BLAS_HH
40#define BLAS_HH
41
42#ifndef F77NAME
43 #define F77NAME(A) A##_
44#endif
45
46#ifndef ATLASNAME
47 #define ATLASNAME(A) ATL_##A
48#endif
49
50namespace blas {
51
52 typedef char character ;
53 typedef int integer ;
54 typedef float single_precision ;
55 typedef double double_precision ;
56
57 #ifdef USE_ATLAS
58
59 /*
60 // # ####### # # #####
61 // # # # # # # # #
62 // # # # # # # #
63 // # # # # # # #####
64 // ####### # # ####### #
65 // # # # # # # # #
66 // # # # ####### # # #####
67 */
68
69 enum ATLAS_ORDER {AtlasRowMajor=101, AtlasColMajor=102 };
70 enum ATLAS_TRANS {AtlasNoTrans=111, AtlasTrans=112,
71 AtlasConjTrans=113, AtlasConj=114};
72 enum ATLAS_UPLO {AtlasUpper=121, AtlasLower=122};
73 enum ATLAS_DIAG {AtlasNonUnit=131, AtlasUnit=132};
74
75 extern "C" {
76 void
77 ATLASNAME(scopy)( const int N,
78 const float *X, const int incX,
79 float *Y, const int incY );
80 void
81 ATLASNAME(dcopy)( const int N,
82 const double *X, const int incX,
83 double *Y, const int incY );
84 void
85 ATLASNAME(sswap)( const int N,
86 float *X, const int incX,
87 float *Y, const int incY );
88 void
89 ATLASNAME(dswap)( const int N,
90 double *X, const int incX,
91 double *Y, const int incY );
92 void
93 ATLASNAME(sscal)( const int N,
94 const float alpha,
95 float *X, const int incX );
96 void
97 ATLASNAME(dscal)( const int N,
98 const double alpha,
99 double *X, const int incX );
100 void
101 ATLASNAME(saxpy)( const int N,
102 const float alpha,
103 const float *X, const int incX,
104 float *Y, const int incY );
105 void
106 ATLASNAME(daxpy)( const int N,
107 const double alpha,
108 const double *X, const int incX,
109 double *Y, const int incY );
110 void
111 ATLASNAME(sger)( const int M, const int N,
112 const float alpha,
113 const float *X, const int incX,
114 const float *Y, const int incY,
115 float *A, const int lda );
116 void
117 ATLASNAME(dger)( const int M, const int N,
118 const double alpha,
119 const double *X, const int incX,
120 const double *Y, const int incY,
121 double *A, const int lda );
122 void
123 ATLASNAME(sgemv)( const enum ATLAS_TRANS TransA,
124 const int M, const int N,
125 const float alpha,
126 const float *A, const int lda,
127 const float *X, const int incX,
128 const float beta,
129 float *Y, const int incY );
130 void
131 ATLASNAME(dgemv)( const enum ATLAS_TRANS TransA,
132 const int M, const int N,
133 const double alpha,
134 const double *A, const int lda,
135 const double *X, const int incX,
136 const double beta,
137 double *Y, const int incY );
138 void
139 ATLASNAME(sgemm)( const enum ATLAS_TRANS TransA,
140 const enum ATLAS_TRANS TransB,
141 const int M, const int N, const int K,
142 const float alpha,
143 const float *A, const int lda,
144 const float *B, const int ldb,
145 const float beta,
146 float *C, const int ldc );
147 void
148 ATLASNAME(dgemm)( const enum ATLAS_TRANS TransA,
149 const enum ATLAS_TRANS TransB,
150 const int M, const int N, const int K,
151 const double alpha,
152 const double *A, const int lda,
153 const double *B, const int ldb,
154 const double beta,
155 double *C, const int ldc );
156 void
157 ATLASNAME(szero)( const int N, float *X, const int incX );
158 void
159 ATLASNAME(dzero)( const int N, double *X, const int incX );
160 float
161 ATLASNAME(snrm2)( const int N, const float *X, const int incX );
162 double
163 ATLASNAME(dnrm2)( const int N, const double *X, const int incX );
164 float
165 ATLASNAME(sasum)( const int N, const float *X, const int incX );
166 double
167 ATLASNAME(dasum)( const int N, const double *X, const int incX );
168 int
169 ATLASNAME(isamax)( const int N, const float *X, const int incX );
170 int
171 ATLASNAME(idamax)( const int N, const double *X, const int incX );
172 float
173 ATLASNAME(sdot)( const int N,
174 const float *X, const int incX,
175 const float *Y, const int incY );
176 double
177 ATLASNAME(ddot)( const int N,
178 const double *X, const int incX,
179 const double *Y, const int incY );
180 }
181
182 inline
183 ATLAS_TRANS
184 trans_atlas(character const TRANS ) {
185 switch ( TRANS ) {
186 case 'N': return AtlasNoTrans ;
187 case 'T': return AtlasTrans ;
188 case 'C': return AtlasConjTrans ;
189 default: throw runtime_error("bad traspose mode selection\n") ;
190 }
191 }
192
193 inline
194 ATLAS_UPLO
195 uplo_atlas(character const UPLO ) {
196 switch ( UPLO ) {
197 case 'U': return AtlasUpper ;
198 case 'L': return AtlasLower ;
199 default: throw runtime_error("bad UP/LOW mode selection\n") ;
200 }
201 }
202
203 inline
204 ATLAS_DIAG
205 diag_atlas(character const DIAG ) {
206 switch ( DIAG ) {
207 case 'U': return AtlasUnit ;
208 case 'N': return AtlasNonUnit ;
209 default: throw runtime_error("bad diagonal mode selection\n") ;
210 }
211 }
212
213 #else
214
215 /*
216 // ##### # ## ####
217 // # # # # # #
218 // ##### # # # ####
219 // # # # ###### #
220 // # # # # # # #
221 // ##### ###### # # ####
222 */
223
224 extern "C" {
225 void
226 F77NAME(scopy)( integer const * N,
227 single_precision const X[],
228 integer const * INCX,
229 single_precision Y[],
230 integer const * INCY ) ;
231 void
232 F77NAME(dcopy)( integer const * N,
233 double_precision const X[],
234 integer const * INCX,
235 double_precision Y[],
236 integer const * INCY ) ;
237 void
238 F77NAME(sswap)( integer const * N,
239 single_precision X[],
240 integer const * INCX,
241 single_precision Y[],
242 integer const * INCY ) ;
243 void
244 F77NAME(dswap)( integer const * N,
245 double_precision X[],
246 integer const * INCX,
247 double_precision Y[],
248 integer const * INCY ) ;
249 void
250 F77NAME(sscal)( integer const * N,
251 single_precision const * S,
252 single_precision X[],
253 integer const * INCX ) ;
254 void
255 F77NAME(dscal)( integer const * N,
256 double_precision const * S,
257 double_precision X[],
258 integer const * INCX ) ;
259 void
260 F77NAME(saxpy)( integer const * N,
261 single_precision const * A,
262 single_precision const X[],
263 integer const * INCX,
264 single_precision Y[],
265 integer const * INCY ) ;
266 void
267 F77NAME(daxpy)( integer const * N,
268 double_precision const * A,
269 double_precision const X[],
270 integer const * INCX,
271 double_precision Y[],
272 integer const * INCY ) ;
273 void
274 F77NAME(strsv)( character const UPLO[],
275 character const TRANS[],
276 character const DIAG[],
277 integer const * N,
278 single_precision const A[],
279 integer const * LDA,
280 single_precision X[],
281 integer const * INCX ) ;
282 void
283 F77NAME(dtrsv)( character const UPLO[],
284 character const TRANS[],
285 character const DIAG[],
286 integer const * N,
287 double_precision const A[],
288 integer const * LDA,
289 double_precision X[],
290 integer const * INCX ) ;
291 void
292 F77NAME(sger)( integer const * M,
293 integer const * N,
294 single_precision const * ALPHA,
295 single_precision const X[],
296 integer const * INCX,
297 single_precision const Y[],
298 integer const * INCY,
299 single_precision A[],
300 integer const * LDA ) ;
301 void
302 F77NAME(dger)( integer const * M,
303 integer const * N,
304 double_precision const * ALPHA,
305 double_precision const X[],
306 integer const * INCX,
307 double_precision const Y[],
308 integer const * INCY,
309 double_precision A[],
310 integer const * LDA ) ;
311 void
312 F77NAME(sgemv)( character const TRANS[],
313 integer const * M,
314 integer const * N,
315 single_precision const * ALPHA,
316 single_precision const A[],
317 integer const * LDA,
318 single_precision const X[],
319 integer const * INCX,
320 single_precision const * BETA,
321 single_precision Y[],
322 integer const * INCY ) ;
323 void
324 F77NAME(dgemv)( character const TRANS[],
325 integer const * M,
326 integer const * N,
327 double_precision const * ALPHA,
328 double_precision const A[],
329 integer const * LDA,
330 double_precision const X[],
331 integer const * INCX,
332 double_precision const * BETA,
333 double_precision Y[],
334 integer const * INCY ) ;
335 void
336 F77NAME(sgemm)( character const TRANSA[],
337 character const TRANSB[],
338 integer const * M,
339 integer const * N,
340 integer const * K,
341 single_precision const * ALPHA,
342 single_precision const A[],
343 integer const * LDA,
344 single_precision const B[],
345 integer const * LDB,
346 single_precision const * BETA,
347 single_precision C[],
348 integer const * LDC ) ;
349 void
350 F77NAME(dgemm)( character const TRANSA[],
351 character const TRANSB[],
352 integer const * M,
353 integer const * N,
354 integer const * K,
355 double_precision const * ALPHA,
356 double_precision const A[],
357 integer const * LDA,
358 double_precision const B[],
359 integer const * LDB,
360 double_precision const * BETA,
361 double_precision C[],
362 integer const * LDC ) ;
363 single_precision
364 F77NAME(snrm2)( integer const * N,
365 single_precision const X[],
366 integer const * INCX ) ;
367
368 double_precision
369 F77NAME(dnrm2)( integer const * N,
370 double_precision const X[],
371 integer const * INCX ) ;
372 single_precision
373 F77NAME(sasum)( integer const * N,
374 single_precision const X[],
375 integer const * INCX ) ;
376 double_precision
377 F77NAME(dasum)( integer const * N,
378 double_precision const X[],
379 integer const * INCX ) ;
380 integer
381 F77NAME(isamax)( integer const * N,
382 single_precision const X[],
383 integer const * INCX ) ;
384 integer
385 F77NAME(idamax)( integer const * N,
386 double_precision const X[],
387 integer const * INCX ) ;
388 single_precision
389 F77NAME(sdot)( integer const * N,
390 single_precision const SX[],
391 integer const * INCX,
392 single_precision const SY[],
393 integer const * INCY ) ;
394
395 double_precision
396 F77NAME(ddot)( integer const * N,
397 double_precision const SX[],
398 integer const * INCX,
399 double_precision const SY[],
400 integer const * INCY ) ;
401 }
402
403 #endif
404
405 /* ____ ____ ___ _ _
406 // | | | |__] \_/
407 // |___ |__| | |
408 */
409 inline
410 void
411 copy( integer const N,
412 single_precision const X[],
413 integer const INCX,
414 single_precision Y[],
415 integer const INCY )
416 #ifdef USE_ATLAS
417 { ATLASNAME(scopy)( N, X, INCX, Y, INCY ) ; }
418 #else
419 { F77NAME(scopy)( &N, X, &INCX, Y, &INCY ) ; }
420 #endif
421
422 inline
423 void
424 copy( integer const N,
425 double_precision const X[],
426 integer const INCX,
427 double_precision Y[],
428 integer const INCY )
429 #ifdef USE_ATLAS
430 { ATLASNAME(dcopy)( N, X, INCX, Y, INCY ) ; }
431 #else
432 { F77NAME(dcopy)( &N, X, &INCX, Y, &INCY ) ; }
433 #endif
434
435 /* __ _
436 // (_ \ / /\ |_)
437 // __) \/\/ /--\ |
438 */
439 inline
440 void
441 swap( integer const N,
442 single_precision X[],
443 integer const INCX,
444 single_precision Y[],
445 integer const INCY )
446 #ifdef USE_ATLAS
447 { ATLASNAME(sswap)( N, X, INCX, Y, INCY ) ; }
448 #else
449 { F77NAME(sswap)( &N, X, &INCX, Y, &INCY ) ; }
450 #endif
451
452 inline
453 void
454 swap( integer const N,
455 double_precision X[],
456 integer const INCX,
457 double_precision Y[],
458 integer const INCY )
459 #ifdef USE_ATLAS
460 { ATLASNAME(dswap)( N, X, INCX, Y, INCY ) ; }
461 #else
462 { F77NAME(dswap)( &N, X, &INCX, Y, &INCY ) ; }
463 #endif
464
465 /* ____ ____ ____ _
466 // [__ | |__| |
467 // ___] |___ | | |___
468 */
469 inline
470 void
471 scal( integer const N,
472 single_precision const S,
473 single_precision X[],
474 integer const INCX )
475 #ifdef USE_ATLAS
476 { ATLASNAME(sscal)(N, S, X, INCX) ; }
477 #else
478 { F77NAME(sscal)(&N, &S, X, &INCX) ; }
479 #endif
480
481 inline
482 void
483 scal( integer const N,
484 double_precision const S,
485 double_precision X[],
486 integer const INCX )
487 #ifdef USE_ATLAS
488 { ATLASNAME(dscal)(N, S, X, INCX) ; }
489 #else
490 { F77NAME(dscal)(&N, &S, X, &INCX) ; }
491 #endif
492
493 /* ____ _ _ ___ _ _
494 // |__| \/ |__] \_/
495 // | | _/\_ | |
496 */
497 inline
498 void
499 axpy( integer const N,
500 single_precision const A,
501 single_precision const X[],
502 integer const INCX,
503 single_precision Y[],
504 integer const INCY )
505 #ifdef USE_ATLAS
506 { ATLASNAME(saxpy)( N, A, X, INCX, Y, INCY ) ; }
507 #else
508 { F77NAME(saxpy)( &N, &A, X, &INCX, Y, &INCY ) ; }
509 #endif
510
511 inline
512 void
513 axpy( integer const N,
514 double_precision const A,
515 double_precision const X[],
516 integer const INCX,
517 double_precision Y[],
518 integer const INCY )
519 #ifdef USE_ATLAS
520 { ATLASNAME(daxpy)( N, A, X, INCX, Y, INCY ) ; }
521 #else
522 { F77NAME(daxpy)( &N, &A, X, &INCX, Y, &INCY ) ; }
523 #endif
524
525 /* __ _ _ _
526 // / |_ |_) / \
527 // /_ |_ | \ \_/
528 */
529 inline
530 void
531 zero( integer const N,
532 single_precision X[],
533 integer const INCX )
534 #ifdef USE_ATLAS
535 { ATLASNAME(szero)( N, X, INCX ) ; }
536 #else
537 { single_precision z = 0 ;
538 integer iz = 0 ;
539 F77NAME(scopy)( &N, &z, &iz, X, &INCX ) ; }
540 #endif
541
542 inline
543 void
544 zero( integer const N,
545 double_precision X[],
546 integer const INCX )
547 #ifdef USE_ATLAS
548 { ATLASNAME(dzero)( N, X, INCX ) ; }
549 #else
550 { double_precision z = 0 ;
551 integer iz = 0 ;
552 F77NAME(dcopy)( &N, &z, &iz, X, &INCX ) ; }
553 #endif
554
555 /* __ _ _
556 // /__ |_ |_)
557 // \_| |_ | \
558 */
559 inline
560 void
561 ger( integer const M,
562 integer const N,
563 single_precision const ALPHA,
564 single_precision const X[],
565 integer const INCX,
566 single_precision const Y[],
567 integer const INCY,
568 single_precision A[],
569 integer const LDA )
570 #ifdef USE_ATLAS
571 { ATLASNAME(sger)( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) ; }
572 #else
573 { F77NAME(sger)( &M, &N, &ALPHA, X, &INCX, Y, &INCY, A, &LDA ) ; }
574 #endif
575
576 inline
577 void
578 ger( integer const M,
579 integer const N,
580 double_precision const ALPHA,
581 double_precision const X[],
582 integer const INCX,
583 double_precision const Y[],
584 integer const INCY,
585 double_precision A[],
586 integer const LDA )
587 #ifdef USE_ATLAS
588 { ATLASNAME(dger)( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) ; }
589 #else
590 { F77NAME(dger)( &M, &N, &ALPHA, X, &INCX, Y, &INCY, A, &LDA ) ; }
591 #endif
592
593 /* __ _
594 // /__ |_ |\/| \ /
595 // \_| |_ | | \/
596 */
597 inline
598 void
599 gemv( character const TRANS[],
600 integer const M,
601 integer const N,
602 single_precision const ALPHA,
603 single_precision const A[],
604 integer const LDA,
605 single_precision const X[],
606 integer const INCX,
607 single_precision const BETA,
608 single_precision Y[],
609 integer const INCY )
610 #ifdef USE_ATLAS
611 { ATLASNAME(sgemv)( trans_atlas(TRANS[0]), M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ; }
612 #else
613 { F77NAME(sgemv)( TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY ) ; }
614 #endif
615
616 inline
617 void
618 gemv( character const TRANS[],
619 integer const M,
620 integer const N,
621 double_precision const ALPHA,
622 double_precision const A[],
623 integer const LDA,
624 double_precision const X[],
625 integer const INCX,
626 double_precision const BETA,
627 double_precision Y[],
628 integer const INCY )
629 #ifdef USE_ATLAS
630 { ATLASNAME(dgemv)( trans_atlas(TRANS[0]), M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ; }
631 #else
632 { F77NAME(dgemv)( TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY ) ; }
633 #endif
634
635 /* ___ _ __
636 // | |_) (_ \ /
637 // | | \ __) \/
638 */
639 inline
640 void
641 trsv( character const UPLO[],
642 character const TRANS[],
643 character const DIAG[],
644 integer const N,
645 single_precision const A[],
646 integer const LDA,
647 single_precision X[],
648 integer const INCX )
649 #ifdef USE_ATLAS
650 { ATLASNAME(strsv)( uplo_atlas(UPLO[0]),
651 trans_atlas(TRANS[0]),
652 diag_atlas(DIAG[0]),
653 N, A, LDA, X, INCX ) ; }
654 #else
655 { F77NAME(strsv)( UPLO, TRANS, DIAG, &N, A, &LDA, X, &INCX ) ; }
656 #endif
657
658 inline
659 void
660 trsv( character const UPLO[],
661 character const TRANS[],
662 character const DIAG[],
663 integer const N,
664 double_precision const A[],
665 integer const LDA,
666 double_precision X[],
667 integer const INCX )
668 #ifdef USE_ATLAS
669 { ATLASNAME(dtrsv)( uplo_atlas(UPLO[0]),
670 trans_atlas(TRANS[0]),
671 diag_atlas(DIAG[0]),
672 N, A, LDA, X, INCX ) ; }
673 #else
674 { F77NAME(dtrsv)( UPLO, TRANS, DIAG, &N, A, &LDA, X, &INCX ) ; }
675 #endif
676
677 /* __ _
678 // /__ |_ |\/| |\/|
679 // \_| |_ | | | |
680 */
681 inline
682 void
683 gemm( character const TRANSA[],
684 character const TRANSB[],
685 integer const M,
686 integer const N,
687 integer const K,
688 single_precision const ALPHA,
689 single_precision const A[],
690 integer const LDA,
691 single_precision const B[],
692 integer const LDB,
693 single_precision const BETA,
694 single_precision C[],
695 integer const LDC )
696 #ifdef USE_ATLAS
697 { ATLASNAME(sgemm)( trans_atlas(TRANSA[0]), trans_atlas(TRANSB[0]),
698 M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) ; }
699 #else
700 { F77NAME(sgemm)( TRANSA, TRANSB,
701 &M, &N, &K,
702 &ALPHA, A, &LDA,
703 B, &LDB,
704 &BETA, C, &LDC ) ; }
705 #endif
706
707 inline
708 void
709 gemm( character const TRANSA[],
710 character const TRANSB[],
711 integer const M,
712 integer const N,
713 integer const K,
714 double_precision const ALPHA,
715 double_precision const A[],
716 integer const LDA,
717 double_precision const B[],
718 integer const LDB,
719 double_precision const BETA,
720 double_precision C[],
721 integer const LDC )
722 #ifdef USE_ATLAS
723 { ATLASNAME(dgemm)( trans_atlas(TRANSA[0]), trans_atlas(TRANSB[0]),
724 M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) ; }
725 #else
726 { F77NAME(dgemm)( TRANSA, TRANSB,
727 &M, &N, &K,
728 &ALPHA, A, &LDA,
729 B, &LDB,
730 &BETA, C, &LDC ) ; }
731 #endif
732
733 /* _ _
734 // |\ | |_) |\/| )
735 // | \| | \ | | /_
736 */
737 inline
738 single_precision
739 nrm2( integer const N,
740 single_precision const X[],
741 integer const INCX )
742 #ifdef USE_ATLAS
743 { return ATLASNAME(snrm2)( N, X, INCX ) ; }
744 #else
745 { return F77NAME(snrm2)( &N, X, &INCX ) ; }
746 #endif
747
748 inline
749 double_precision
750 nrm2( integer const N,
751 double_precision const X[],
752 integer const INCX )
753 #ifdef USE_ATLAS
754 { return ATLASNAME(dnrm2)( N, X, INCX ) ; }
755 #else
756 { return F77NAME(dnrm2)( &N, X, &INCX ) ; }
757 #endif
758
759 /* __
760 // /\ (_ | ||\/|
761 // /--\__)|_|| |
762 */
763 inline
764 single_precision
765 asum(integer const N,
766 single_precision const X[],
767 integer const INCX)
768 #ifdef USE_ATLAS
769 { return ATLASNAME(sasum)( N, X, INCX ) ; }
770 #else
771 { return F77NAME(sasum)( &N, X, &INCX ) ; }
772 #endif
773
774 inline
775 double_precision
776 asum(integer const N,
777 double_precision const X[],
778 integer const INCX)
779 #ifdef USE_ATLAS
780 { return ATLASNAME(dasum)( N, X, INCX ) ; }
781 #else
782 { return F77NAME(dasum)( &N, X, &INCX ) ; }
783 #endif
784
785 /*
786 // /\ |\/| /\ \/
787 // /--\ | | /--\ /\
788 */
789 inline
790 single_precision
791 amax(integer const N,
792 single_precision const X[],
793 integer const INCX)
794 #ifdef USE_ATLAS
795 { return X[ATLASNAME(isamax)( N, X, INCX )-1] ; }
796 #else
797 { return X[F77NAME(isamax)( &N, X, &INCX )-1] ; }
798 #endif
799
800 inline
801 double_precision
802 amax(integer const N,
803 double_precision const X[],
804 integer const INCX)
805 #ifdef USE_ATLAS
806 { return X[ATLASNAME(idamax)( N, X, INCX )-1] ; }
807 #else
808 { return X[F77NAME(idamax)( &N, X, &INCX )-1] ; }
809 #endif
810
811 /* _ _ ___
812 // | \ / \ |
813 // |_/ \_/ |
814 */
815 inline
816 single_precision
817 dot(integer const N,
818 single_precision const SX[],
819 integer const INCX,
820 single_precision const SY[],
821 integer const INCY)
822 #ifdef USE_ATLAS
823 { return ATLASNAME(sdot)( N, SX, INCX, SY, INCY ) ; }
824 #else
825 { return F77NAME(sdot)( &N, SX, &INCX, SY, &INCY ) ; }
826 #endif
827
828 inline
829 double_precision
830 dot(integer const N,
831 double_precision const SX[],
832 integer const INCX,
833 double_precision const SY[],
834 integer const INCY)
835 #ifdef USE_ATLAS
836 { return ATLASNAME(ddot)( N, SX, INCX, SY, INCY ) ; }
837 #else
838 { return F77NAME(ddot)( &N, SX, &INCX, SY, &INCY ) ; }
839 #endif
840
841} // end namespace blas
842
843#endif
844
845///
846/// eof: Blas.hh
847///
Lapack C/FORTRAN interaface
File bvpbetter/Lapack.hh
:open:
1/*--------------------------------------------------------------------------*\
2 | |
3 | This program is free software; you can redistribute it and/or modify |
4 | it under the terms of the GNU General Public License as published by |
5 | the Free Software Foundation; either version 2, or (at your option) |
6 | any later version. |
7 | |
8 | This program is distributed in the hope that it will be useful, |
9 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
10 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
11 | GNU General Public License for more details. |
12 | |
13 | You should have received a copy of the GNU General Public License |
14 | along with this program; if not, write to the Free Software |
15 | Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
16 | |
17 | Copyright (C) 2008 |
18 | |
19 | , __ , __ |
20 | /|/ \ /|/ \ |
21 | | __/ _ ,_ | __/ _ ,_ |
22 | | \|/ / | | | | \|/ / | | | |
23 | |(__/|__/ |_/ \_/|/|(__/|__/ |_/ \_/|/ |
24 | /| /| |
25 | \| \| |
26 | |
27 | Enrico Bertolazzi |
28 | Dipartimento di Ingegneria Meccanica e Strutturale |
29 | Universita` degli Studi di Trento |
30 | Via Mesiano 77, I-38050 Trento, Italy |
31 | email: enrico.bertolazzi@unitn.it |
32 | |
33\*--------------------------------------------------------------------------*/
34
35///
36/// file: Lapack.hh
37///
38
39#ifndef LAPACK_HH
40#define LAPACK_HH
41
42#ifndef F77NAME
43 #define F77NAME(A) A##_
44#endif
45
46#ifndef ATLASNAME
47 #define ATLASNAME(A) ATL_##A
48#endif
49
50namespace lapack {
51
52 typedef char character ;
53 typedef int integer ;
54 typedef float single_precision ;
55 typedef double double_precision ;
56
57 #ifdef USE_ATLAS
58
59 /*
60 // # ####### # # #####
61 // # # # # # # # #
62 // # # # # # # #
63 // # # # # # # #####
64 // ####### # # ####### #
65 // # # # # # # # #
66 // # # # ####### # # #####
67 */
68
69 enum ATLAS_ORDER {AtlasRowMajor=101, AtlasColMajor=102 };
70 enum ATLAS_TRANS {AtlasNoTrans=111, AtlasTrans=112,
71 AtlasConjTrans=113, AtlasConj=114};
72 enum ATLAS_UPLO {AtlasUpper=121, AtlasLower=122};
73 enum ATLAS_DIAG {AtlasNonUnit=131, AtlasUnit=132};
74
75 extern "C" {
76 void
77 ATLASNAME(sgetrs)( const enum ATLAS_ORDER Order,
78 const enum ATLAS_TRANS Trans,
79 const int N,
80 const int NRHS,
81 const float *A, const int lda,
82 const int *ipiv,
83 float *B, const int ldb );
84 void
85 ATLASNAME(dgetrs)( const enum ATLAS_ORDER Order,
86 const enum ATLAS_TRANS Trans,
87 const int N,
88 const int NRHS,
89 const double *A, const int lda,
90 const int *ipiv,
91 double *B, const int ldb );
92 int
93 ATLASNAME(sgetrf)( const enum ATLAS_ORDER Order,
94 const int M, const int N,
95 float *A, const int lda,
96 int *ipiv );
97 int
98 ATLASNAME(dgetrf)( const enum ATLAS_ORDER Order,
99 const int M, const int N,
100 double *A, const int lda,
101 int *ipiv );
102 void
103 ATLASNAME(strsv)( const enum ATLAS_UPLO Uplo,
104 const enum ATLAS_TRANS TransA,
105 const enum ATLAS_DIAG Diag,
106 const int N,
107 const float *A, const int lda,
108 float *X, const int incX );
109 void
110 ATLASNAME(dtrsv)( const enum ATLAS_UPLO Uplo,
111 const enum ATLAS_TRANS TransA,
112 const enum ATLAS_DIAG Diag,
113 const int N,
114 const double *A, const int lda,
115 double *X, const int incX );
116 void
117 ATLASNAME(szero)( const int N, float *X, const int incX );
118 void
119 ATLASNAME(dzero)( const int N, double *X, const int incX );
120
121 void
122 ATLASNAME(sgezero)( const int M, const int N, float *C, const int ldc );
123 void
124 ATLASNAME(dgezero)( const int M, const int N, double *C, const int ldc );
125
126 void
127 ATLASNAME(sgecopy)( const int M, const int N,
128 const float *A, const int lda,
129 float *C, const int ldc );
130 void
131 ATLASNAME(dgecopy)( const int M, const int N,
132 const double *A, const int lda,
133 double *C, const int ldc );
134 }
135
136 inline
137 ATLAS_TRANS
138 trans_atlas(character const TRANS ) {
139 switch ( TRANS ) {
140 case 'N': return AtlasNoTrans ;
141 case 'T': return AtlasTrans ;
142 case 'C': return AtlasConjTrans ;
143 default: throw runtime_error("bad traspose mode selection\n") ;
144 }
145 }
146
147 inline
148 ATLAS_UPLO
149 uplo_atlas(character const UPLO ) {
150 switch ( UPLO ) {
151 case 'U': return AtlasUpper ;
152 case 'L': return AtlasLower ;
153 default: throw runtime_error("bad UP/LOW mode selection\n") ;
154 }
155 }
156
157 inline
158 ATLAS_DIAG
159 diag_atlas(character const DIAG ) {
160 switch ( DIAG ) {
161 case 'U': return AtlasUnit ;
162 case 'N': return AtlasNonUnit ;
163 default: throw runtime_error("bad diagonal mode selection\n") ;
164 }
165 }
166
167 #else
168
169 /*
170 // # ## ##### ## #### # #
171 // # # # # # # # # # # #
172 // # # # # # # # # ####
173 // # ###### ##### ###### # # #
174 // # # # # # # # # # #
175 // ###### # # # # # #### # #
176 */
177
178 extern "C" {
179 void
180 F77NAME(dgetrf)( integer const * N,
181 integer const * M,
182 double_precision A[],
183 integer const * LDA,
184 integer IPIV[],
185 integer * INFO ) ;
186 void
187 F77NAME(sgetrs)( character const TRANS[],
188 integer const * N,
189 integer const * NRHS,
190 single_precision const A[],
191 integer const * LDA,
192 integer const IPIV[],
193 single_precision B[],
194 integer const * LDB,
195 integer * INFO ) ;
196 void
197 F77NAME(dgetrs)( character const TRANS[],
198 integer const * N,
199 integer const * NRHS,
200 double_precision const A[],
201 integer const * LDA,
202 integer const IPIV[],
203 double_precision B[],
204 integer const * LDB,
205 integer * INFO ) ;
206 void
207 F77NAME(sgetrf)( integer const * N,
208 integer const * M,
209 single_precision A[],
210 integer const * LDA,
211 integer IPIV[],
212 integer * INFO ) ;
213 void
214 F77NAME(dgetrf)( integer const * N,
215 integer const * M,
216 double_precision A[],
217 integer const * LDA,
218 integer IPIV[],
219 integer * INFO ) ;
220 void
221 F77NAME(sgetrs)( character const TRANS[],
222 integer const * N,
223 integer const * NRHS,
224 single_precision const A[],
225 integer const * LDA,
226 integer const IPIV[],
227 single_precision B[],
228 integer const * LDB,
229 integer * INFO ) ;
230 void
231 F77NAME(dgetrs)( character const TRANS[],
232 integer const * N,
233 integer const * NRHS,
234 double_precision const A[],
235 integer const * LDA,
236 integer const IPIV[],
237 double_precision B[],
238 integer const * LDB,
239 integer * INFO ) ;
240 void
241 F77NAME(sgesv)( integer const * N,
242 integer const * NRHS,
243 single_precision A[],
244 integer const * LDA,
245 integer IPIV[],
246 single_precision B[],
247 integer const * LDB,
248 integer * INFO ) ;
249 void
250 F77NAME(dgesv)( integer const * N,
251 integer const * NRHS,
252 double_precision A[],
253 integer const * LDA,
254 integer IPIV[],
255 double_precision B[],
256 integer const * LDB,
257 integer * INFO ) ;
258 void
259 F77NAME(strsv)( character const UPLO[],
260 character const TRANS[],
261 character const DIAG[],
262 integer const * N,
263 single_precision const A[],
264 integer const * LDA,
265 single_precision X[],
266 integer const * INCX ) ;
267 void
268 F77NAME(dtrsv)( character const UPLO[],
269 character const TRANS[],
270 character const DIAG[],
271 integer const * N,
272 double_precision const A[],
273 integer const * LDA,
274 double_precision X[],
275 integer const * INCX ) ;
276 void
277 F77NAME(slaset)( character const UPLO[],
278 integer const * M,
279 integer const * N,
280 single_precision const * ALPHA,
281 single_precision const * BETA,
282 single_precision A[],
283 integer const * LDA ) ;
284 void
285 F77NAME(dlaset)( character const UPLO[],
286 integer const * M,
287 integer const * N,
288 double_precision const * ALPHA,
289 double_precision const * BETA,
290 double_precision A[],
291 integer const * LDA ) ;
292 void
293 F77NAME(slacpy)( character const UPLO[],
294 integer const * M,
295 integer const * N,
296 single_precision const A[],
297 integer const * LDA,
298 single_precision B[],
299 integer const * LDB) ;
300 void
301 F77NAME(dlacpy)( character const UPLO[],
302 integer const * M,
303 integer const * N,
304 double_precision const A[],
305 integer const * LDA,
306 double_precision B[],
307 integer const * LDB ) ;
308 }
309
310 #endif
311
312 /* ____ ____ ___ _ _
313 // | | | |__] \_/
314 // |___ |__| | |
315 */
316 inline
317 void
318 gecopy( integer const M,
319 integer const N,
320 single_precision const A[],
321 integer const LDA,
322 single_precision B[],
323 integer const LDB )
324 #ifdef USE_ATLAS
325 { ATLASNAME(sgecopy)( M, N, A, LDA, B, LDB ) ; }
326 #else
327 { F77NAME(slacpy)( "A", &M, &N, A, &LDA, B, &LDB ) ; }
328 #endif
329
330 inline
331 void
332 gecopy( integer const M,
333 integer const N,
334 double_precision const A[],
335 integer const LDA,
336 double_precision B[],
337 integer const LDB )
338 #ifdef USE_ATLAS
339 { ATLASNAME(dgecopy)( M, N, A, LDA, B, LDB ) ; }
340 #else
341 { F77NAME(dlacpy)( "A", &M, &N, A, &LDA, B, &LDB ) ; }
342 #endif
343
344 /* __ _ _ _
345 // / |_ |_) / \
346 // /_ |_ | \ \_/
347 */
348 inline
349 void
350 gezero( integer const M,
351 integer const N,
352 single_precision A[],
353 integer const LDA )
354 #ifdef USE_ATLAS
355 { ATLASNAME(sgezero)( M, N, A, LDA ) ; }
356 #else
357 { single_precision zero = 0 ;
358 F77NAME(slaset)( "A", &M, &N, &zero, &zero, A, &LDA ) ; }
359 #endif
360
361 inline
362 void
363 gezero( integer const M,
364 integer const N,
365 double_precision A[],
366 integer const LDA )
367 #ifdef USE_ATLAS
368 { ATLASNAME(dgezero)( M, N, A, LDA ) ; }
369 #else
370 { double_precision zero = 0 ;
371 F77NAME(dlaset)( "A", &M, &N, &zero, &zero, A, &LDA ) ; }
372 #endif
373
374 /* ___ _ __
375 // | |_) (_ \ /
376 // | | \ __) \/
377 */
378 inline
379 void
380 trsv( character const UPLO[],
381 character const TRANS[],
382 character const DIAG[],
383 integer const N,
384 single_precision const A[],
385 integer const LDA,
386 single_precision X[],
387 integer const INCX )
388 #ifdef USE_ATLAS
389 { ATLASNAME(strsv)( uplo_atlas(UPLO[0]),
390 trans_atlas(TRANS[0]),
391 diag_atlas(DIAG[0]),
392 N, A, LDA, X, INCX ) ; }
393 #else
394 { F77NAME(strsv)( UPLO, TRANS, DIAG, &N, A, &LDA, X, &INCX ) ; }
395 #endif
396
397 inline
398 void
399 trsv( character const UPLO[],
400 character const TRANS[],
401 character const DIAG[],
402 integer const N,
403 double_precision const A[],
404 integer const LDA,
405 double_precision X[],
406 integer const INCX )
407 #ifdef USE_ATLAS
408 { ATLASNAME(dtrsv)( uplo_atlas(UPLO[0]),
409 trans_atlas(TRANS[0]),
410 diag_atlas(DIAG[0]),
411 N, A, LDA, X, INCX ) ; }
412 #else
413 { F77NAME(dtrsv)( UPLO, TRANS, DIAG, &N, A, &LDA, X, &INCX ) ; }
414 #endif
415
416 /* __ _ ___ _ _
417 // /__ |_ | |_) |_
418 // \_| |_ | | \ |
419 */
420 inline
421 integer
422 getrf( integer const N,
423 integer const M,
424 single_precision A[],
425 integer const LDA,
426 integer IPIV[] )
427 #ifdef USE_ATLAS
428 { return ATLASNAME(sgetrf)( AtlasColMajor, N, M, A, LDA, IPIV ) ; }
429 #else
430 { integer INFO ;
431 F77NAME(sgetrf)( &N, &M, A, &LDA, IPIV, &INFO ) ;
432 return INFO ; }
433 #endif
434
435 inline
436 integer
437 getrf( integer const N,
438 integer const M,
439 double_precision A[],
440 integer const LDA,
441 integer IPIV[] )
442 #ifdef USE_ATLAS
443 { return ATLASNAME(dgetrf)( AtlasColMajor, N, M, A, LDA, IPIV ) ; }
444 #else
445 { integer INFO ;
446 F77NAME(dgetrf)( &N, &M, A, &LDA, IPIV, &INFO ) ;
447 return INFO ; }
448 #endif
449
450 /* __ _ ___ _ __
451 // /__ |_ | |_) (_
452 // \_| |_ | | \ __)
453 */
454 inline
455 integer
456 getrs( character const TRANS[],
457 integer const N,
458 integer const NRHS,
459 single_precision const A[],
460 integer const LDA,
461 integer const IPIV[],
462 single_precision B[],
463 integer const LDB )
464 #ifdef USE_ATLAS
465 { ATLASNAME(sgetrs)( AtlasColMajor, trans_atlas(TRANS[0]), N, NRHS, A, LDA, IPIV, B, LDB ) ;
466 return 0 ; }
467 #else
468 { integer INFO ;
469 F77NAME(sgetrs)( TRANS, &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO ) ;
470 return INFO ; }
471 #endif
472
473 inline
474 integer
475 getrs( character const TRANS[],
476 integer const N,
477 integer const NRHS,
478 double_precision const A[],
479 integer const LDA,
480 integer const IPIV[],
481 double_precision B[],
482 integer const LDB )
483 #ifdef USE_ATLAS
484 { ATLASNAME(dgetrs)( AtlasColMajor, trans_atlas(TRANS[0]), N, NRHS, A, LDA, IPIV, B, LDB ) ;
485 return 0 ; }
486 #else
487 { integer INFO ;
488 F77NAME(dgetrs)( TRANS, &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO ) ;
489 return INFO ; }
490 #endif
491
492 /* __ _ __
493 // /__ |_ (_ \ /
494 // \_| |_ __) \/
495 */
496 inline
497 integer
498 gesv( integer const N,
499 integer const NRHS,
500 single_precision A[],
501 integer const LDA,
502 integer IPIV[],
503 single_precision B[],
504 integer const LDB )
505 #ifdef USE_ATLAS
506 { integer INFO = ATLASNAME(sgetrf)(AtlasColMajor, N, N, A, LDA, IPIV);
507 if ( INFO == 0 ) ATLASNAME(sgetrs)(AtlasColMajor, AtlasNoTrans,
508 N, NRHS, A, LDA, IPIV, B, LDB) ;
509 return INFO ; }
510 #else
511 { integer INFO ;
512 F77NAME(sgesv)( &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO ) ;
513 return INFO ; }
514 #endif
515
516 inline
517 integer
518 gesv( integer const N,
519 integer const NRHS,
520 double_precision A[],
521 integer const LDA,
522 integer IPIV[],
523 double_precision B[],
524 integer const LDB )
525 #ifdef USE_ATLAS
526 { integer INFO = ATLASNAME(dgetrf)(AtlasColMajor, N, N, A, LDA, IPIV);
527 if ( INFO == 0 ) ATLASNAME(dgetrs)(AtlasColMajor, AtlasNoTrans,
528 N, NRHS, A, LDA, IPIV, B, LDB) ;
529 return INFO ; }
530 #else
531 { integer INFO ;
532 F77NAME(dgesv)( &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO ) ;
533 return INFO ; }
534 #endif
535
536 extern "C" {
537
538 void
539 F77NAME(sgttrf)( integer const * N,
540 single_precision DL[],
541 single_precision D[],
542 single_precision DU[],
543 single_precision DU2[],
544 integer IPIV[],
545 integer * INFO ) ;
546
547 void
548 F77NAME(dgttrf)( integer const * N,
549 double_precision DL[],
550 double_precision D[],
551 double_precision DU[],
552 double_precision DU2[],
553 integer IPIV[],
554 integer * INFO ) ;
555
556 void
557 F77NAME(sgttrs)( character TRANS[],
558 integer const * N,
559 integer const * NRHS,
560 single_precision const DL[],
561 single_precision const D[],
562 single_precision const DU[],
563 single_precision const DU2[],
564 integer const IPIV[],
565 single_precision B[],
566 integer const * LDB,
567 integer * INFO ) ;
568
569 void
570 F77NAME(dgttrs)( character TRANS[],
571 integer const * N,
572 integer const * NRHS,
573 double_precision const DL[],
574 double_precision const D[],
575 double_precision const DU[],
576 double_precision const DU2[],
577 integer const IPIV[],
578 double_precision B[],
579 integer const * LDB,
580 integer * INFO ) ;
581
582 void
583 F77NAME(sgtsv)( integer const * N,
584 integer const * NRHS,
585 single_precision DL[],
586 single_precision D[],
587 single_precision DU[],
588 single_precision B[],
589 integer const * LDB,
590 integer * INFO ) ;
591
592 void
593 F77NAME(dgtsv)( integer const * N,
594 integer const * NRHS,
595 double_precision DL[],
596 double_precision D[],
597 double_precision DU[],
598 double_precision B[],
599 integer const * LDB,
600 integer * INFO ) ;
601 }
602
603 /* __ ___ ___ _ _
604 // /__ | | |_) |_
605 // \_| | | | \ |
606 */
607
608 inline
609 integer
610 gttrf( integer const N,
611 single_precision DL[],
612 single_precision D[],
613 single_precision DU[],
614 single_precision DU2[],
615 integer IPIV[] )
616 { integer INFO ;
617 F77NAME(sgttrf)( &N, DL, D, DU, DU2, IPIV, &INFO) ;
618 return INFO ; }
619
620 inline
621 integer
622 gttrf( integer const N,
623 double_precision DL[],
624 double_precision D[],
625 double_precision DU[],
626 double_precision DU2[],
627 integer IPIV[] )
628 { integer INFO ;
629 F77NAME(dgttrf)( &N, DL, D, DU, DU2, IPIV, &INFO) ;
630 return INFO ; }
631
632 /* __ ___ ___ _ __
633 // /__ | | |_) (_
634 // \_| | | | \ __)
635 */
636
637 inline
638 integer
639 gttrs( character TRANS[],
640 integer const N,
641 integer const NRHS,
642 single_precision const DL[],
643 single_precision const D[],
644 single_precision const DU[],
645 single_precision const DU2[],
646 integer const IPIV[],
647 single_precision B[],
648 integer const LDB )
649 { integer INFO ;
650 F77NAME(sgttrs)( TRANS, &N, &NRHS, DL, D, DU, DU2, IPIV, B, &LDB, &INFO) ;
651 return INFO ; }
652
653 inline
654 integer
655 gttrs( character TRANS[],
656 integer const N,
657 integer const NRHS,
658 double_precision const DL[],
659 double_precision const D[],
660 double_precision const DU[],
661 double_precision const DU2[],
662 integer const IPIV[],
663 double_precision B[],
664 integer const LDB )
665 { integer INFO ;
666 F77NAME(dgttrs)( TRANS, &N, &NRHS, DL, D, DU, DU2, IPIV, B, &LDB, &INFO) ;
667 return INFO ; }
668
669 /* __ ___ __
670 // /__ | (_ \ /
671 // \_| | __) \/
672 */
673
674 inline
675 integer
676 gtsv( integer const N,
677 integer const NRHS,
678 single_precision DL[],
679 single_precision D[],
680 single_precision DU[],
681 single_precision B[],
682 integer const LDB )
683 { integer INFO ;
684 F77NAME(sgtsv)( &N, &NRHS, DL, D, DU, B, &LDB, &INFO ) ;
685 return INFO ; }
686
687 inline
688 integer
689 gtsv( integer const N,
690 integer const NRHS,
691 double_precision DL[],
692 double_precision D[],
693 double_precision DU[],
694 double_precision B[],
695 integer const LDB )
696 { integer INFO ;
697 F77NAME(dgtsv)( &N, &NRHS, DL, D, DU, B, &LDB, &INFO ) ;
698 return INFO ; }
699
700} // end namespace lapack
701
702#endif
703
704///
705/// eof: Lapack.hh
706///
The lesson in a zip file
The zip file lesson7.zip