MAGMA  magma-1.4.0
Matrix Algebra on GPU and Multicore Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
testing_zgetrf.cpp File Reference
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cuda_runtime_api.h>
#include <cublas.h>
#include "flops.h"
#include "magma.h"
#include "magma_lapack.h"
#include "testings.h"
Include dependency graph for testing_zgetrf.cpp:

Go to the source code of this file.

Functions

void init_matrix (int m, int n, magmaDoubleComplex *h_A, magma_int_t lda)
 
double get_residual (magma_int_t m, magma_int_t n, magmaDoubleComplex *A, magma_int_t lda, magma_int_t *ipiv)
 
double get_LU_error (magma_int_t M, magma_int_t N, magmaDoubleComplex *LU, magma_int_t lda, magma_int_t *ipiv)
 
int main (int argc, char **argv)
 

Function Documentation

double get_LU_error ( magma_int_t  M,
magma_int_t  N,
magmaDoubleComplex *  LU,
magma_int_t  lda,
magma_int_t ipiv 
)

Definition at line 102 of file testing_zgetrf.cpp.

References A, blasf77_zgemm, init_matrix(), L, lapackf77_zlacpy, lapackf77_zlange, lapackf77_zlaswp, MAGMA_Z_MAKE, MAGMA_Z_ONE, MAGMA_Z_SUB, MAGMA_Z_ZERO, MagmaLowerStr, MagmaUpperStr, min, TESTING_FREE, and TESTING_MALLOC.

105 {
106  magma_int_t min_mn = min(M,N);
107  magma_int_t ione = 1;
108  magma_int_t i, j;
109  magmaDoubleComplex alpha = MAGMA_Z_ONE;
110  magmaDoubleComplex beta = MAGMA_Z_ZERO;
111  magmaDoubleComplex *A, *L, *U;
112  double work[1], matnorm, residual;
113 
114  TESTING_MALLOC( A, magmaDoubleComplex, lda*N );
115  TESTING_MALLOC( L, magmaDoubleComplex, M*min_mn );
116  TESTING_MALLOC( U, magmaDoubleComplex, min_mn*N );
117  memset( L, 0, M*min_mn*sizeof(magmaDoubleComplex) );
118  memset( U, 0, min_mn*N*sizeof(magmaDoubleComplex) );
119 
120  // set to original A
121  init_matrix( M, N, A, lda );
122  lapackf77_zlaswp( &N, A, &lda, &ione, &min_mn, ipiv, &ione);
123 
124  // copy LU to L and U, and set diagonal to 1
125  lapackf77_zlacpy( MagmaLowerStr, &M, &min_mn, LU, &lda, L, &M );
126  lapackf77_zlacpy( MagmaUpperStr, &min_mn, &N, LU, &lda, U, &min_mn );
127  for(j=0; j<min_mn; j++)
128  L[j+j*M] = MAGMA_Z_MAKE( 1., 0. );
129 
130  matnorm = lapackf77_zlange("f", &M, &N, A, &lda, work);
131 
132  blasf77_zgemm("N", "N", &M, &N, &min_mn,
133  &alpha, L, &M, U, &min_mn, &beta, LU, &lda);
134 
135  for( j = 0; j < N; j++ ) {
136  for( i = 0; i < M; i++ ) {
137  LU[i+j*lda] = MAGMA_Z_SUB( LU[i+j*lda], A[i+j*lda] );
138  }
139  }
140  residual = lapackf77_zlange("f", &M, &N, LU, &lda, work);
141 
142  TESTING_FREE(A);
143  TESTING_FREE(L);
144  TESTING_FREE(U);
145 
146  return residual / (matnorm * N);
147 }
#define min(a, b)
Definition: common_magma.h:86
#define MAGMA_Z_MAKE(r, i)
Definition: magma.h:123
#define MAGMA_Z_SUB(a, b)
Definition: magma.h:127
void init_matrix(int m, int n, magmaFloatComplex *h_A, magma_int_t lda)
int magma_int_t
Definition: magmablas.h:12
#define TESTING_MALLOC(__ptr, __type, __size)
Definition: testings.h:34
#define TESTING_FREE(__ptr)
Definition: testings.h:54
#define lapackf77_zlacpy
Definition: magma_zlapack.h:73
#define L(i, j)
#define lapackf77_zlaswp
Definition: magma_zlapack.h:85
#define LU(i, j)
#define MAGMA_Z_ZERO
Definition: magma.h:131
#define MagmaLowerStr
Definition: magma.h:85
#define A(i, j)
Definition: cprint.cpp:16
#define blasf77_zgemm
Definition: magma_zlapack.h:33
#define lapackf77_zlange
Definition: magma_zlapack.h:75
#define MAGMA_Z_ONE
Definition: magma.h:132
#define MagmaUpperStr
Definition: magma.h:84

Here is the call graph for this function:

double get_residual ( magma_int_t  m,
magma_int_t  n,
magmaDoubleComplex *  A,
magma_int_t  lda,
magma_int_t ipiv 
)

Definition at line 43 of file testing_zgetrf.cpp.

References blasf77_zcopy, blasf77_zgemv, init_matrix(), lapackf77_zgetrs, lapackf77_zlange, lapackf77_zlarnv, magma_strerror(), MAGMA_Z_NEG_ONE, MAGMA_Z_ONE, TESTING_FREE, and TESTING_MALLOC.

47 {
48  if ( m != n ) {
49  printf( "\nERROR: residual check defined only for square matrices\n" );
50  return -1;
51  }
52 
53  const magmaDoubleComplex c_one = MAGMA_Z_ONE;
54  const magmaDoubleComplex c_neg_one = MAGMA_Z_NEG_ONE;
55  const magma_int_t ione = 1;
56 
57  // this seed should be DIFFERENT than used in init_matrix
58  // (else x is column of A, so residual can be exactly zero)
59  magma_int_t ISEED[4] = {0,0,0,2};
60  magma_int_t info = 0;
61  magmaDoubleComplex *x, *b;
62 
63  // initialize RHS
64  TESTING_MALLOC( x, magmaDoubleComplex, n );
65  TESTING_MALLOC( b, magmaDoubleComplex, n );
66  lapackf77_zlarnv( &ione, ISEED, &n, b );
67  blasf77_zcopy( &n, b, &ione, x, &ione );
68 
69  // solve Ax = b
70  lapackf77_zgetrs( "Notrans", &n, &ione, A, &lda, ipiv, x, &n, &info );
71  if (info != 0)
72  printf("lapackf77_zgetrs returned error %d: %s.\n",
73  (int) info, magma_strerror( info ));
74 
75  // reset to original A
76  init_matrix( m, n, A, lda );
77 
78  // compute r = Ax - b, saved in b
79  blasf77_zgemv( "Notrans", &m, &n, &c_one, A, &lda, x, &ione, &c_neg_one, b, &ione );
80 
81  // compute residual |Ax - b| / (n*|A|*|x|)
82  double norm_x, norm_A, norm_r, work[1];
83  norm_A = lapackf77_zlange( "F", &m, &n, A, &lda, work );
84  norm_r = lapackf77_zlange( "F", &n, &ione, b, &n, work );
85  norm_x = lapackf77_zlange( "F", &n, &ione, x, &n, work );
86 
87  //printf( "r=\n" ); magma_zprint( 1, n, b, 1 );
88 
89  TESTING_FREE( x );
90  TESTING_FREE( b );
91 
92  //printf( "r=%.2e, A=%.2e, x=%.2e, n=%d\n", norm_r, norm_A, norm_x, n );
93  return norm_r / (n * norm_A * norm_x);
94 }
#define MAGMA_Z_NEG_ONE
Definition: magma.h:134
void init_matrix(int m, int n, magmaFloatComplex *h_A, magma_int_t lda)
int magma_int_t
Definition: magmablas.h:12
#define TESTING_MALLOC(__ptr, __type, __size)
Definition: testings.h:34
#define TESTING_FREE(__ptr)
Definition: testings.h:54
#define A(i, j)
Definition: cprint.cpp:16
#define blasf77_zcopy
Definition: magma_zlapack.h:24
const char * magma_strerror(magma_err_t error)
Definition: error.cpp:227
#define lapackf77_zlange
Definition: magma_zlapack.h:75
#define MAGMA_Z_ONE
Definition: magma.h:132
#define blasf77_zgemv
Definition: magma_zlapack.h:34
#define lapackf77_zgetrs
Definition: magma_zlapack.h:65
#define lapackf77_zlarnv
Definition: magma_zlapack.h:81

Here is the call graph for this function:

void init_matrix ( int  m,
int  n,
magmaDoubleComplex *  h_A,
magma_int_t  lda 
)

Definition at line 29 of file testing_zgetrf.cpp.

References lapackf77_zlarnv.

30 {
31  magma_int_t ione = 1;
32  magma_int_t ISEED[4] = {0,0,0,1};
33  magma_int_t n2 = lda*n;
34  lapackf77_zlarnv( &ione, ISEED, &n2, h_A );
35 }
int magma_int_t
Definition: magmablas.h:12
#define h_A(i, j)
#define lapackf77_zlarnv
Definition: magma_zlapack.h:81
int main ( int  argc,
char **  argv 
)

Definition at line 153 of file testing_zgetrf.cpp.

References magma_opts::check, FLOPS_ZGETRF, get_LU_error(), get_residual(), h_A, init_matrix(), magma_opts::lapack, lapackf77_dlamch, lapackf77_zgetrf, magma_strerror(), magma_wtime(), magma_zgetrf(), min, magma_opts::msize, magma_opts::ngpu, magma_opts::niter, magma_opts::nsize, magma_opts::ntest, parse_opts(), TESTING_FINALIZE, TESTING_FREE, TESTING_HOSTALLOC, TESTING_HOSTFREE, TESTING_INIT, TESTING_MALLOC, and magma_opts::tolerance.

154 {
155  TESTING_INIT();
156 
157  real_Double_t gflops, gpu_perf, gpu_time, cpu_perf=0, cpu_time=0;
158  double error;
159  magmaDoubleComplex *h_A;
160  magma_int_t *ipiv;
161  magma_int_t M, N, n2, lda, ldda, info, min_mn;
162  magma_int_t status = 0;
163 
164  magma_opts opts;
165  parse_opts( argc, argv, &opts );
166 
167  double tol = opts.tolerance * lapackf77_dlamch("E");
168 
169  printf("ngpu %d\n", (int) opts.ngpu );
170  if ( opts.check == 2 ) {
171  printf(" M N CPU GFlop/s (sec) GPU GFlop/s (sec) |Ax-b|/(N*|A|*|x|)\n");
172  }
173  else {
174  printf(" M N CPU GFlop/s (sec) GPU GFlop/s (sec) |PA-LU|/(N*|A|)\n");
175  }
176  printf("=========================================================================\n");
177  for( int i = 0; i < opts.ntest; ++i ) {
178  for( int iter = 0; iter < opts.niter; ++iter ) {
179  M = opts.msize[i];
180  N = opts.nsize[i];
181  min_mn = min(M, N);
182  lda = M;
183  n2 = lda*N;
184  ldda = ((M+31)/32)*32;
185  gflops = FLOPS_ZGETRF( M, N ) / 1e9;
186 
187  TESTING_MALLOC( ipiv, magma_int_t, min_mn );
188  TESTING_HOSTALLOC( h_A, magmaDoubleComplex, n2 );
189 
190  /* =====================================================================
191  Performs operation using LAPACK
192  =================================================================== */
193  if ( opts.lapack ) {
194  init_matrix( M, N, h_A, lda );
195 
196  cpu_time = magma_wtime();
197  lapackf77_zgetrf(&M, &N, h_A, &lda, ipiv, &info);
198  cpu_time = magma_wtime() - cpu_time;
199  cpu_perf = gflops / cpu_time;
200  if (info != 0)
201  printf("lapackf77_zgetrf returned error %d: %s.\n",
202  (int) info, magma_strerror( info ));
203  }
204 
205  /* ====================================================================
206  Performs operation using MAGMA
207  =================================================================== */
208  init_matrix( M, N, h_A, lda );
209 
210  gpu_time = magma_wtime();
211  magma_zgetrf( M, N, h_A, lda, ipiv, &info);
212  gpu_time = magma_wtime() - gpu_time;
213  gpu_perf = gflops / gpu_time;
214  if (info != 0)
215  printf("magma_zgetrf returned error %d: %s.\n",
216  (int) info, magma_strerror( info ));
217 
218  /* =====================================================================
219  Check the factorization
220  =================================================================== */
221  if ( opts.lapack ) {
222  printf("%5d %5d %7.2f (%7.2f) %7.2f (%7.2f)",
223  (int) M, (int) N, cpu_perf, cpu_time, gpu_perf, gpu_time );
224  }
225  else {
226  printf("%5d %5d --- ( --- ) %7.2f (%7.2f)",
227  (int) M, (int) N, gpu_perf, gpu_time );
228  }
229  if ( opts.check == 2 ) {
230  error = get_residual( M, N, h_A, lda, ipiv );
231  printf(" %8.2e%s\n", error, (error < tol ? "" : " failed"));
232  status |= ! (error < tol);
233  }
234  else if ( opts.check ) {
235  error = get_LU_error( M, N, h_A, lda, ipiv );
236  printf(" %8.2e%s\n", error, (error < tol ? "" : " failed"));
237  status |= ! (error < tol);
238  }
239  else {
240  printf(" --- \n");
241  }
242 
243  TESTING_FREE( ipiv );
244  TESTING_HOSTFREE( h_A );
245  }
246  if ( opts.niter > 1 ) {
247  printf( "\n" );
248  }
249  }
250 
252  return status;
253 }
void parse_opts(int argc, char **argv, magma_opts *opts)
#define min(a, b)
Definition: common_magma.h:86
#define TESTING_HOSTFREE(__ptr)
Definition: testings.h:57
magma_int_t ntest
Definition: testings.h:124
#define lapackf77_zgetrf
Definition: magma_zlapack.h:64
void init_matrix(int m, int n, magmaFloatComplex *h_A, magma_int_t lda)
magma_int_t niter
Definition: testings.h:138
#define TESTING_INIT()
Definition: testings.h:19
int magma_int_t
Definition: magmablas.h:12
#define TESTING_MALLOC(__ptr, __type, __size)
Definition: testings.h:34
#define TESTING_FREE(__ptr)
Definition: testings.h:54
double get_LU_error(magma_int_t M, magma_int_t N, cuDoubleComplex *A, magma_int_t lda, cuDoubleComplex *LU, magma_int_t *IPIV)
#define h_A(i, j)
float get_residual(magma_int_t m, magma_int_t n, magmaFloatComplex *A, magma_int_t lda, magma_int_t *ipiv)
#define lapackf77_dlamch
Definition: magma_lapack.h:27
magma_int_t magma_zgetrf(magma_int_t m, magma_int_t n, cuDoubleComplex *A, magma_int_t lda, magma_int_t *ipiv, magma_int_t *info)
magma_int_t ngpu
Definition: testings.h:137
int check
Definition: testings.h:147
magma_int_t ldda
double magma_wtime(void)
Definition: timer.cpp:110
magma_int_t nsize[MAX_NTEST]
Definition: testings.h:126
const char * magma_strerror(magma_err_t error)
Definition: error.cpp:227
#define TESTING_HOSTALLOC(__ptr, __type, __size)
Definition: testings.h:41
#define TESTING_FINALIZE()
Definition: testings.h:29
double real_Double_t
Definition: magma_types.h:27
#define FLOPS_ZGETRF(__m, __n)
Definition: flops.h:223
magma_int_t msize[MAX_NTEST]
Definition: testings.h:125
double tolerance
Definition: testings.h:144
int lapack
Definition: testings.h:148

Here is the call graph for this function: