test_check_correctness_template.F90 18.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), formerly known as
!      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!
!    More information can be found here:
!    http://elpa.mpcdf.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA 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 Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
! Author: A. Marek, MPCDF

44
    function check_correctness_evp_numeric_residuals_&
45
46
47
    &MATH_DATATYPE&
    &_&
    &PRECISION&
Pavel Kus's avatar
Pavel Kus committed
48
    & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol) result(status)
49
      implicit none
50
#include "../../src/general/precision_kinds.F90"
51
      integer(kind=ik)                 :: status
Pavel Kus's avatar
Pavel Kus committed
52
      integer(kind=ik), intent(in)     :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
53
54
55
56
      MATH_DATATYPE(kind=rck), intent(in)     :: as(:,:), z(:,:)
      real(kind=rk)                 :: ev(:)
      MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
      MATH_DATATYPE(kind=rck)                :: xc
57

58
#ifndef WITH_MPI
59
60
#if REALCASE == 1
      real(kind=rck)                   :: dnrm2, snrm2
61
62
#endif
#if COMPLEXCASE == 1
63
      complex(kind=rck)                :: zdotc, cdotc
64
#endif /* COMPLEXCASE */
65
#endif
66
67
68

      integer(kind=ik)                 :: sc_desc(:)

Pavel Kus's avatar
Pavel Kus committed
69
      integer(kind=ik)                 :: i, rowLocal, colLocal
70
      real(kind=rck)                   :: err, errmax
71
72
73

      integer :: mpierr

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
      ! tolerance for the residual test for different math type/precision setups
      real(kind=rk), parameter       :: tol_res_real_double      = 5e-12_rk
      real(kind=rk), parameter       :: tol_res_real_single      = 3e-3_rk
      real(kind=rk), parameter       :: tol_res_complex_double   = 5e-12_rk
      real(kind=rk), parameter       :: tol_res_complex_single   = 3e-3_rk
      real(kind=rk), parameter       :: tol_res                  = tol_res_&
                                                                          &MATH_DATATYPE&
                                                                          &_&
                                                                          &PRECISION
      ! tolerance for the orthogonality test for different math type/precision setups
      real(kind=rk), parameter       :: tol_orth_real_double     = 5e-12_rk
      real(kind=rk), parameter       :: tol_orth_real_single     = 9e-4_rk
      real(kind=rk), parameter       :: tol_orth_complex_double  = 5e-12_rk
      real(kind=rk), parameter       :: tol_orth_complex_single  = 9e-4_rk
      real(kind=rk), parameter       :: tol_orth                 = tol_orth_&
                                                                          &MATH_DATATYPE&
                                                                          &_&
                                                                          &PRECISION

93
94
95
96
97
98
99
      status = 0

      ! 1. Residual (maximum of || A*Zi - Zi*EVi ||)
      ! tmp1 =  A * Z
      ! as is original stored matrix, Z are the EVs

#ifdef WITH_MPI
100
101
      call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, &
                  z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
102
#else /* WITH_MPI */
103
      call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na)
104
105
106
107
108
109
#endif /* WITH_MPI */


      ! tmp2 = Zi*EVi
      tmp2(:,:) = z(:,:)
      do i=1,nev
Pavel Kus's avatar
Pavel Kus committed
110
        xc = ev(i)
111
#ifdef WITH_MPI
Pavel Kus's avatar
Pavel Kus committed
112
113
114
        call p&
            &BLAS_CHAR&
            &scal(na, xc, tmp2, 1, i, sc_desc, 1)
115
#else /* WITH_MPI */
Pavel Kus's avatar
Pavel Kus committed
116
117
        call BLAS_CHAR&
            &scal(na,xc,tmp2(:,i),1)
118
119
120
121
122
123
124
#endif /* WITH_MPI */
      enddo

      !  tmp1 = A*Zi - Zi*EVi
      tmp1(:,:) =  tmp1(:,:) - tmp2(:,:)

      ! Get maximum norm of columns of tmp1
125
      errmax = 0.0_rk
126
127
128

      do i=1,nev
#if REALCASE == 1
129
        err = 0.0_rk
130
#ifdef WITH_MPI
Pavel Kus's avatar
Pavel Kus committed
131
        call scal_PRECISION_NRM2(na, err, tmp1, 1, i, sc_desc, 1)
132
#else /* WITH_MPI */
Pavel Kus's avatar
Pavel Kus committed
133
        err = PRECISION_NRM2(na,tmp1(1,i),1)
134
135
136
137
138
139
140
#endif /* WITH_MPI */
        errmax = max(errmax, err)
#endif /* REALCASE */

#if COMPLEXCASE == 1
        xc = 0
#ifdef WITH_MPI
Pavel Kus's avatar
Pavel Kus committed
141
        call scal_PRECISION_DOTC(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
142
#else /* WITH_MPI */
Pavel Kus's avatar
Pavel Kus committed
143
        xc = PRECISION_DOTC(na,tmp1,1,tmp1,1)
144
#endif /* WITH_MPI */
145
        errmax = max(errmax, sqrt(real(xc,kind=REAL_DATATYPE)))
146
147
148
149
150
151
#endif /* COMPLEXCASE */
      enddo

      ! Get maximum error norm over all processors
      err = errmax
#ifdef WITH_MPI
152
      call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
153
154
155
#else /* WITH_MPI */
      errmax = err
#endif /* WITH_MPI */
156
      if (myid==0) print *,'Results of numerical residual checks:'
157
      if (myid==0) print *,'Error Residual     :',errmax
158
      if (nev .ge. 2) then
159
        if (errmax .gt. tol_res .or. errmax .eq. 0.0_rk) then
160
161
162
          status = 1
        endif
      else
163
        if (errmax .gt. tol_res) then
164
165
          status = 1
        endif
166
167
168
169
170
171
172
      endif

      ! 2. Eigenvector orthogonality

      ! tmp1 = Z**T * Z
      tmp1 = 0
#ifdef WITH_MPI
173
174
      call scal_PRECISION_GEMM(BLAS_TRANS_OR_CONJ, 'N', nev, nev, na, ONE, z, 1, 1, &
                        sc_desc, z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
175
#else /* WITH_MPI */
176
      call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,z,na,ZERO,tmp1,na)
177
#endif /* WITH_MPI */
Pavel Kus's avatar
Pavel Kus committed
178
179
      !TODO for the C interface, not all information is passed (zeros instead)
      !TODO than this part of the test cannot be done
180
      !TODO either we will not have this part of test at all, or it will be improved
Pavel Kus's avatar
Pavel Kus committed
181
182
      if(nblk > 0) then
        ! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
183
        err = 0.0_rk
Pavel Kus's avatar
Pavel Kus committed
184
185
        do i=1, nev
          if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
186
             err = max(err, abs(tmp1(rowLocal,colLocal) - 1.0_rk))
Pavel Kus's avatar
Pavel Kus committed
187
188
189
190
191
192
193
194
195
           endif
        end do
#ifdef WITH_MPI
        call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else /* WITH_MPI */
        errmax = err
#endif /* WITH_MPI */
        if (myid==0) print *,'Maximal error in eigenvector lengths:',errmax
      end if
196

Pavel Kus's avatar
Pavel Kus committed
197
      ! Second, find the maximal error in the whole Z**T * Z matrix (its diference from identity matrix)
198
199
200
      ! Initialize tmp2 to unit matrix
      tmp2 = 0
#ifdef WITH_MPI
Pavel Kus's avatar
Pavel Kus committed
201
      call scal_PRECISION_LASET('A', nev, nev, ZERO, ONE, tmp2, 1, 1, sc_desc)
202
#else /* WITH_MPI */
Pavel Kus's avatar
Pavel Kus committed
203
      call PRECISION_LASET('A',nev,nev,ZERO,ONE,tmp2,na)
204
205
206
207
208
209
210
211
#endif /* WITH_MPI */

      !      ! tmp1 = Z**T * Z - Unit Matrix
      tmp1(:,:) =  tmp1(:,:) - tmp2(:,:)

      ! Get maximum error (max abs value in tmp1)
      err = maxval(abs(tmp1))
#ifdef WITH_MPI
212
      call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
213
214
215
216
#else /* WITH_MPI */
      errmax = err
#endif /* WITH_MPI */
      if (myid==0) print *,'Error Orthogonality:',errmax
217
      if (nev .ge. 2) then
218
        if (errmax .gt. tol_orth .or. errmax .eq. 0.0_rk) then
219
220
221
          status = 1
        endif
      else
222
        if (errmax .gt. tol_orth) then
223
224
225
          status = 1
        endif
      endif
226
227
228
229
230
231
    end function


#if REALCASE == 1

#ifdef DOUBLE_PRECISION_REAL
232
    !c> int check_correctness_evp_numeric_residuals_real_double_f(int na, int nev, int na_rows, int na_cols,
233
234
235
    !c>                                         double *as, double *z, double *ev,
    !c>                                         int sc_desc[9], int myid);
#else
236
    !c> int check_correctness_evp_numeric_residuals_real_single_f(int na, int nev, int na_rows, int na_cols,
237
238
239
240
241
242
243
244
    !c>                                         float *as, float *z, float *ev,
    !c>                                         int sc_desc[9], int myid);
#endif

#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
245
    !c> int check_correctness_evp_numeric_residuals_complex_double_f(int na, int nev, int na_rows, int na_cols,
246
247
248
    !c>                                         complex double *as, complex double *z, double *ev,
    !c>                                         int sc_desc[9], int myid);
#else
249
    !c> int check_correctness_evp_numeric_residuals_complex_single_f(int na, int nev, int na_rows, int na_cols,
Andreas Marek's avatar
Andreas Marek committed
250
    !c>                                         complex float *as, complex float *z, float *ev,
251
252
253
254
    !c>                                         int sc_desc[9], int myid);
#endif
#endif /* COMPLEXCASE */

255
function check_correctness_evp_numeric_residuals_&
256
257
258
259
&MATH_DATATYPE&
&_&
&PRECISION&
&_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid) result(status) &
260
      bind(C,name="check_correctness_evp_numeric_residuals_&
261
262
263
264
265
266
267
268
      &MATH_DATATYPE&
      &_&
      &PRECISION&
      &_f")

      use iso_c_binding

      implicit none
269
#include "../../src/general/precision_kinds.F90"
270
271
272

      integer(kind=c_int)            :: status
      integer(kind=c_int), value     :: na, nev, myid, na_rows, na_cols
273
274
      MATH_DATATYPE(kind=rck)     :: as(1:na_rows,1:na_cols), z(1:na_rows,1:na_cols)
      real(kind=rck)    :: ev(1:na)
275
276
      integer(kind=c_int)            :: sc_desc(1:9)

Pavel Kus's avatar
Pavel Kus committed
277
278
279
      ! TODO: I did not want to add all the variables to the C interface as well
      ! TODO: I think that we should find a better way to pass this information
      ! TODO: to all the functions anyway (get it from sc_desc, pass elpa_t, etc..)
280
      status = check_correctness_evp_numeric_residuals_&
281
282
283
      &MATH_DATATYPE&
      &_&
      &PRECISION&
Pavel Kus's avatar
Pavel Kus committed
284
      & (na, nev, as, z, ev, sc_desc, 0, myid, 0, 0, 0, 0)
285
286
287

    end function

Andreas Marek's avatar
Andreas Marek committed
288
289
290
291
292
293
294
    function check_correctness_eigenvalues_toeplitz_&
    &MATH_DATATYPE&
    &_&
    &PRECISION&
    & (na, diagonalElement, subdiagonalElement, ev, z, myid) result(status)
      use iso_c_binding
      implicit none
295
#include "../../src/general/precision_kinds.F90"
Andreas Marek's avatar
Andreas Marek committed
296
297
298

      integer               :: status, ii, j, myid
      integer, intent(in)   :: na
299
300
301
      real(kind=rck) :: diagonalElement, subdiagonalElement
      real(kind=rck) :: ev_analytic(na), ev(na)
      MATH_DATATYPE(kind=rck) :: z(:,:)
Andreas Marek's avatar
Andreas Marek committed
302
303

#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
304
      real(kind=rck), parameter   :: pi = 3.141592653589793238462643383279_c_double
Andreas Marek's avatar
Andreas Marek committed
305
#else
306
      real(kind=rck), parameter   :: pi = 3.1415926535897932_c_float
Andreas Marek's avatar
Andreas Marek committed
307
#endif
308
      real(kind=rck)              :: tmp, maxerr
309
      integer                     :: loctmp
Andreas Marek's avatar
Andreas Marek committed
310
311
312
313
      status = 0

     ! analytic solution
     do ii=1, na
314
315
316
       ev_analytic(ii) = diagonalElement + 2.0_rk * &
                         subdiagonalElement *cos( pi*real(ii,kind=rk)/ &
                         real(na+1,kind=rk) )
Andreas Marek's avatar
Andreas Marek committed
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
     enddo

     ! sort analytic solution:

     ! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
     ! a proper sorting algorithmus might be implemented here

     tmp    = minval(ev_analytic)
     loctmp = minloc(ev_analytic, 1)

     ev_analytic(loctmp) = ev_analytic(1)
     ev_analytic(1) = tmp
     do ii=2, na
       tmp = ev_analytic(ii)
       do j= ii, na
         if (ev_analytic(j) .lt. tmp) then
           tmp    = ev_analytic(j)
           loctmp = j
         endif
       enddo
       ev_analytic(loctmp) = ev_analytic(ii)
       ev_analytic(ii) = tmp
     enddo

     ! compute a simple error max of eigenvalues
     maxerr = 0.0
     maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)

#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
346
     if (maxerr .gt. 8.e-13_c_double .or. maxerr .eq. 0.0_c_double) then
Andreas Marek's avatar
Andreas Marek committed
347
#else
348
     if (maxerr .gt. 8.e-4_c_float .or. maxerr .eq. 0.0_c_float) then
Andreas Marek's avatar
Andreas Marek committed
349
350
351
352
353
354
355
#endif
       status = 1
       if (myid .eq. 0) then
         print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
       endif
     endif
    end function
356

357
358
359
360
361
362
363
364
365
366
    function check_correctness_cholesky_&
    &MATH_DATATYPE&
    &_&
    &PRECISION&
    & (na, a, as, na_rows, sc_desc, myid) result(status)
      implicit none
#include "../../src/general/precision_kinds.F90"
      integer(kind=ik)                 :: status
      integer(kind=ik), intent(in)     :: na, myid, na_rows

367
368
369
      MATH_DATATYPE(kind=rck), intent(in)       :: a(:,:), as(:,:)
      MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
      real(kind=rk)                   :: norm, normmax
370
371

#ifdef WITH_MPI
372
373
374
      real(kind=rck)                   :: p&
                                           &BLAS_CHAR&
                                           &lange
375
#else /* WITH_MPI */
376
377
      real(kind=rck)                   :: BLAS_CHAR&
                                          &lange
378
379
380
381
382
383
384
#endif /* WITH_MPI */

      integer(kind=ik)                 :: sc_desc(:)
      real(kind=rck)                   :: err, errmax
      integer :: mpierr

      status = 0
385
      tmp1(:,:) = 0.0_rck
386
387
388
389
390


#if REALCASE == 1
      ! tmp1 = a**T
#ifdef WITH_MPI
391
392
393
      call p&
          &BLAS_CHAR&
          &tran(na, na, 1.0_rck, a, 1, 1, sc_desc, 0.0_rck, tmp1, 1, 1, sc_desc)
394
395
396
397
398
399
400
401
#else /* WITH_MPI */
      tmp1 = transpose(a)
#endif /* WITH_MPI */
#endif /* REALCASE == 1 */

#if COMPLEXCASE == 1
      ! tmp1 = a**H
#ifdef WITH_MPI
402
403
404
      call p&
            &BLAS_CHAR&
            &tranc(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
405
406
407
#else /* WITH_MPI */
      tmp1 = transpose(conjg(a))
#endif /* WITH_MPI */
408
#endif /* COMPLEXCASE == 1 */
409
410
411

      ! tmp2 = a * a**T
#ifdef WITH_MPI
412
413
414
415
      call p&
            &BLAS_CHAR&
            &gemm("N","N", na, na, na, ONE, a, 1, 1, sc_desc, tmp1, 1, 1, &
               sc_desc, ZERO, tmp2, 1, 1, sc_desc)
416
#else /* WITH_MPI */
417
418
      call BLAS_CHAR&
                    &gemm("N","N", na, na, na, ONE, a, na, tmp1, na, ZERO, tmp2, na)
419
420
421
422
#endif /* WITH_MPI */

      ! compare tmp2 with original matrix
      tmp2(:,:) = tmp2(:,:) - as(:,:)
423

424
#ifdef WITH_MPI
425
426
427
      norm = p&
              &BLAS_CHAR&
              &lange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
428
#else /* WITH_MPI */
429
430
      norm = BLAS_CHAR&
             &lange("M", na, na, tmp2, na_rows, tmp1)
431
432
433
434
#endif /* WITH_MPI */


#ifdef WITH_MPI
435
      call mpi_allreduce(norm,normmax,1,MPI_REAL_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr)
436
437
438
439
440
441
442
443
#else /* WITH_MPI */
      normmax = norm
#endif /* WITH_MPI */

      if (myid .eq. 0) then
        print *," Maximum error of result: ", normmax
      endif

444
#if REALCASE == 1
445
#ifdef DOUBLE_PRECISION_REAL
446
!      if (normmax .gt. 5e-12_rk8 .or. normmax .eq. 0.0_rk8) then
447
448
449
450
      if (normmax .gt. 5e-12_rk8) then
        status = 1
      endif
#else
451
452
!      if (normmax .gt. 5e-4_rk4 .or. normmax .eq. 0.0_rk4) then
      if (normmax .gt. 5e-4_rk4 ) then
453
454
455
        status = 1
      endif
#endif
456
#endif
457

458
459
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
460
461
!      if (normmax .gt. 5e-11_rk8 .or. normmax .eq. 0.0_rk8) then
      if (normmax .gt. 5e-11_rk8 ) then
462
463
464
        status = 1
      endif
#else
465
!      if (normmax .gt. 5e-3_rk4 .or. normmax .eq. 0.0_rk4) then
466
467
468
469
470
      if (normmax .gt. 5e-3_rk4) then
        status = 1
      endif
#endif
#endif
471
472
    end function

473
474
475
476
477
478
479
480
481
    function check_correctness_hermitian_multiply_&
    &MATH_DATATYPE&
    &_&
    &PRECISION&
    & (na, a, b, c, na_rows, sc_desc, myid) result(status)
      implicit none
#include "../../src/general/precision_kinds.F90"
      integer(kind=ik)                 :: status
      integer(kind=ik), intent(in)     :: na, myid, na_rows
482
483
      MATH_DATATYPE(kind=rck), intent(in)       :: a(:,:), b(:,:), c(:,:)
      MATH_DATATYPE(kind=rck), dimension(size(a,dim=1),size(a,dim=2)) :: tmp1, tmp2
484
485
486
      real(kind=rck)                   :: norm, normmax

#ifdef WITH_MPI
487
488
489
      real(kind=rck)                   :: p&
                                           &BLAS_CHAR&
                                           &lange
490
#else /* WITH_MPI */
491
492
      real(kind=rck)                   :: BLAS_CHAR&
                                          &lange
493
494
495
496
497
498
499
#endif /* WITH_MPI */

      integer(kind=ik)                 :: sc_desc(:)
      real(kind=rck)                   :: err, errmax
      integer :: mpierr

      status = 0
500
      tmp1(:,:) = ZERO
501
502
503
504

#if REALCASE == 1
      ! tmp1 = a**T
#ifdef WITH_MPI
505
506
507
      call p&
            &BLAS_CHAR&
            &tran(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
508
509
510
511
512
513
514
515
516
#else /* WITH_MPI */
      tmp1 = transpose(a)
#endif /* WITH_MPI */

#endif /* REALCASE == 1 */

#if COMPLEXCASE == 1
      ! tmp1 = a**H
#ifdef WITH_MPI
517
518
519
      call p&
            &BLAS_CHAR&
            &tranc(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
520
521
522
#else /* WITH_MPI */
      tmp1 = transpose(conjg(a))
#endif /* WITH_MPI */
523
#endif /* COMPLEXCASE == 1 */
524
525
526

   ! tmp2 = tmp1 * b
#ifdef WITH_MPI
527
528
529
530
   call p&
         &BLAS_CHAR&
         &gemm("N","N", na, na, na, ONE, tmp1, 1, 1, sc_desc, b, 1, 1, &
               sc_desc, ZERO, tmp2, 1, 1, sc_desc)
531
#else
532
533
   call BLAS_CHAR&
        &gemm("N","N", na, na, na, ONE, tmp1, na, b, na, ZERO, tmp2, na)
534
535
536
537
538
539
#endif

      ! compare tmp2 with c
      tmp2(:,:) = tmp2(:,:) - c(:,:)

#ifdef WITH_MPI
540
541
542
      norm = p&
              &BLAS_CHAR&
              &lange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
543
#else /* WITH_MPI */
544
545
      norm = BLAS_CHAR&
             &lange("M", na, na, tmp2, na_rows, tmp1)
546
547
548
#endif /* WITH_MPI */

#ifdef WITH_MPI
549
      call mpi_allreduce(norm,normmax,1,MPI_REAL_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr)
550
551
552
553
554
555
556
557
558
#else /* WITH_MPI */
      normmax = norm
#endif /* WITH_MPI */

      if (myid .eq. 0) then
        print *," Maximum error of result: ", normmax
      endif

#ifdef DOUBLE_PRECISION_REAL
559
      if (normmax .gt. 5e-11_rk8 ) then
560
561
562
        status = 1
      endif
#else
563
      if (normmax .gt. 5e-3_rk4 ) then
564
565
566
567
        status = 1
      endif
#endif

568
569
570
571
572
573
574
575
576
#ifdef DOUBLE_PRECISION_COMPLEX
      if (normmax .gt. 5e-11_rk8 ) then
        status = 1
      endif
#else
      if (normmax .gt. 5e-3_rk4 ) then
        status = 1
      endif
#endif
577
578
579
    end function


580
581


582
! vim: syntax=fortran