compute template splitted, tridiag file changed to resemble more the real case

Conflicts:
	src/elpa1_compute_complex_template.X90

Makefile.am

@@ -886,6 +886,9 @@ EXTRA_DIST = \
   src/elpa1_trans_ev_real_template.X90 \
   src/elpa1_tridiag_real_template.X90 \
   src/elpa1_compute_complex_template.X90 \
+  src/elpa1_tools_complex_template.X90 \
+  src/elpa1_trans_ev_complex_template.X90 \
+  src/elpa1_tridiag_complex_template.X90 \
   src/elpa2_compute_real_template.X90 \
   src/elpa2_compute_complex_template.X90 \
   src/precision_macros.h \

src/elpa1_tools_complex_template.X90

+    subroutine hh_transform_complex_PRECISION(alpha, xnorm_sq, xf, tau)
+
+      ! Similar to LAPACK routine ZLARFP, but uses ||x||**2 instead of x(:)
+      ! and returns the factor xf by which x has to be scaled.
+      ! It also hasn't the special handling for numbers < 1.d-300 or > 1.d150
+      ! since this would be expensive for the parallel implementation.
+      use precision
+#ifdef HAVE_DETAILED_TIMINGS
+      use timings
+#else
+      use timings_dummy
+#endif
+      implicit none
+      complex(kind=COMPLEX_DATATYPE), intent(inout) :: alpha
+      real(kind=REAL_DATATYPE), intent(in)       :: xnorm_sq
+      complex(kind=COMPLEX_DATATYPE), intent(out)   :: xf, tau
+
+      real(kind=REAL_DATATYPE)                   :: ALPHR, ALPHI, BETA
+      
+      call timer%start("hh_transform_complex" // PRECISION_SUFFIX)
+
+      ALPHR = real( ALPHA, kind=REAL_DATATYPE )
+      ALPHI = PRECISION_IMAG( ALPHA )
+      if ( XNORM_SQ==0. .AND. ALPHI==0. ) then
+
+        if ( ALPHR>=0. ) then
+          TAU = 0.
+        else
+          TAU = 2.
+          ALPHA = -ALPHA
+        endif
+        XF = 0.
+
+      else
+
+        BETA = SIGN( SQRT( ALPHR**2 + ALPHI**2 + XNORM_SQ ), ALPHR )
+        ALPHA = ALPHA + BETA
+        IF ( BETA<0 ) THEN
+          BETA = -BETA
+          TAU = -ALPHA / BETA
+        ELSE
+          ALPHR = ALPHI * (ALPHI/real( ALPHA , kind=KIND_PRECISION))
+          ALPHR = ALPHR + XNORM_SQ/real( ALPHA, kind=KIND_PRECISION )
+
+          TAU = PRECISION_CMPLX( ALPHR/BETA, -ALPHI/BETA )
+          ALPHA = PRECISION_CMPLX( -ALPHR, ALPHI )
+        END IF
+        XF = CONST_REAL_1_0/ALPHA
+        ALPHA = BETA
+      endif
+
+      call timer%stop("hh_transform_complex" // PRECISION_SUFFIX)
+
+    end subroutine hh_transform_complex_PRECISION

src/elpa1_trans_ev_complex_template.X90

+#if 0
+!    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
+!
+!    This particular source code file contains additions, changes and
+!    enhancements authored by Intel Corporation which is not part of
+!    the ELPA consortium.
+!
+!    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.
+!
+!
+! ELPA1 -- Faster replacements for ScaLAPACK symmetric eigenvalue routines
+!
+! Copyright of the original code rests with the authors inside the ELPA
+! consortium. The copyright of any additional modifications shall rest
+! with their original authors, but shall adhere to the licensing terms
+! distributed along with the original code in the file "COPYING".
+#endif
+
+
+!> \brief Transforms the eigenvectors of a tridiagonal matrix back
+!>                     to the eigenvectors of the original matrix
+!>                     (like Scalapack Routine PDORMTR)
+!>
+!  Parameters
+!
+!> \param na          Order of matrix a, number of rows of matrix q
+!>
+!> \param nqc         Number of columns of matrix q
+!>
+!> \param a(lda,matrixCols)  Matrix containing the Householder vectors (i.e. matrix a after tridiag_real)
+!>                           Distribution is like in Scalapack.
+!>
+!> \param lda         Leading dimension of a
+!>
+!> \param tau(na)     Factors of the Householder vectors
+!>
+!> \param q           On input: Eigenvectors of tridiagonal matrix
+!>                    On output: Transformed eigenvectors
+!>                    Distribution is like in Scalapack.
+!>    
+!> \param ldq         Leading dimension of q
+!>
+!> \param nblk        blocksize of cyclic distribution, must be the same in both directions!
+!>
+!> \param matrixCols  local columns of matrix a and q
+!>
+!> \param mpi_comm_rows        MPI-Communicator for rows
+!>
+!> \param mpi_comm_cols        MPI-Communicator for columns
+!>
+!> \param useGPU      If true,  GPU version of the subroutine will be used
+!>
+    subroutine trans_ev_complex_PRECISION(na, nqc, a, lda, tau, q, ldq, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols)
+
+#ifdef HAVE_DETAILED_TIMINGS
+      use timings
+#else
+      use timings_dummy
+#endif
+      use precision
+      implicit none
+
+      integer(kind=ik)              ::  na, nqc, lda, ldq, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols
+      complex(kind=COMPLEX_DATATYPE)              ::  tau(na)
+#ifdef USE_ASSUMED_SIZE
+      complex(kind=COMPLEX_DATATYPE)              :: a(lda,*), q(ldq,*)
+#else
+      complex(kind=COMPLEX_DATATYPE)              ::  a(lda,matrixCols), q(ldq,matrixCols)
+#endif
+      integer(kind=ik)              :: max_stored_rows
+#ifdef DOUBLE_PRECISION_COMPLEX
+      complex(kind=ck8), parameter   :: CZERO = (0.0_rk8,0.0_rk8), CONE = (1.0_rk8,0.0_rk8)
+#else
+      complex(kind=ck4), parameter   :: CZERO = (0.0_rk4,0.0_rk4), CONE = (1.0_rk4,0.0_rk4)
+#endif
+      integer(kind=ik)              :: my_prow, my_pcol, np_rows, np_cols, mpierr
+      integer(kind=ik)              :: totalblocks, max_blocks_row, max_blocks_col, max_local_rows, max_local_cols
+      integer(kind=ik)              :: l_cols, l_rows, l_colh, nstor
+      integer(kind=ik)              :: istep, i, n, nc, ic, ics, ice, nb, cur_pcol
+
+      complex(kind=COMPLEX_DATATYPE), allocatable :: tmp1(:), tmp2(:), hvb(:), hvm(:,:)
+      complex(kind=COMPLEX_DATATYPE), allocatable :: tmat(:,:), h1(:), h2(:)
+      integer(kind=ik)              :: istat
+      character(200)                :: errorMessage
+
+      call timer%start("trans_ev_complex" // PRECISION_SUFFIX)
+#ifdef HAVE_DETAILED_TIMINGS
+      call timer%start("mpi_communication")
+#endif
+
+      call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
+      call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)
+      call mpi_comm_rank(mpi_comm_cols,my_pcol,mpierr)
+      call mpi_comm_size(mpi_comm_cols,np_cols,mpierr)
+#ifdef HAVE_DETAILED_TIMINGS
+      call timer%stop("mpi_communication")
+#endif
+
+      totalblocks = (na-1)/nblk + 1
+      max_blocks_row = (totalblocks-1)/np_rows + 1
+      max_blocks_col = ((nqc-1)/nblk)/np_cols + 1  ! Columns of q!
+
+      max_local_rows = max_blocks_row*nblk
+      max_local_cols = max_blocks_col*nblk
+
+      max_stored_rows = (63/nblk+1)*nblk
+
+      allocate(tmat(max_stored_rows,max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating tmat "//errorMessage
+       stop
+      endif
+
+      allocate(h1(max_stored_rows*max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating h1 "//errorMessage
+       stop
+      endif
+
+      allocate(h2(max_stored_rows*max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating h2 "//errorMessage
+       stop
+      endif
+
+      allocate(tmp1(max_local_cols*max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating tmp1 "//errorMessage
+       stop
+      endif
+
+      allocate(tmp2(max_local_cols*max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating tmp2 "//errorMessage
+       stop
+      endif
+
+      allocate(hvb(max_local_rows*nblk), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating hvb "//errorMessage
+       stop
+      endif
+
+      allocate(hvm(max_local_rows,max_stored_rows), stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when allocating hvm "//errorMessage
+       stop
+      endif
+
+      hvm = 0   ! Must be set to 0 !!!
+      hvb = 0   ! Safety only
+
+      l_cols = local_index(nqc, my_pcol, np_cols, nblk, -1) ! Local columns of q
+
+      nstor = 0
+
+      ! In the complex case tau(2) /= 0
+      if (my_prow == prow(1, nblk, np_rows)) then
+        q(1,1:l_cols) = q(1,1:l_cols)*(CONE-tau(2))
+      endif
+
+      do istep=1,na,nblk
+
+        ics = MAX(istep,3)
+        ice = MIN(istep+nblk-1,na)
+        if (ice<ics) cycle
+
+        cur_pcol = pcol(istep, nblk, np_cols)
+
+        nb = 0
+        do ic=ics,ice
+
+          l_colh = local_index(ic  , my_pcol, np_cols, nblk, -1) ! Column of Householder vector
+          l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector
+
+
+          if (my_pcol==cur_pcol) then
+            hvb(nb+1:nb+l_rows) = a(1:l_rows,l_colh)
+            if (my_prow==prow(ic-1, nblk, np_rows)) then
+              hvb(nb+l_rows) = 1.
+            endif
+          endif
+
+          nb = nb+l_rows
+        enddo
+
+#ifdef WITH_MPI
+#ifdef HAVE_DETAILED_TIMINGS
+        call timer%start("mpi_communication")
+#endif
+
+        if (nb>0) &
+           call MPI_Bcast(hvb, nb, MPI_COMPLEX_PRECISION, cur_pcol, mpi_comm_cols, mpierr)
+#ifdef HAVE_DETAILED_TIMINGS
+        call timer%stop("mpi_communication")
+#endif
+
+#endif /* WITH_MPI */
+        nb = 0
+        do ic=ics,ice
+          l_rows = local_index(ic-1, my_prow, np_rows, nblk, -1) ! # rows of Householder vector
+          hvm(1:l_rows,nstor+1) = hvb(nb+1:nb+l_rows)
+          nstor = nstor+1
+          nb = nb+l_rows
+        enddo
+
+        ! Please note: for smaller matix sizes (na/np_rows<=256), a value of 32 for nstor is enough!
+        if (nstor+nblk>max_stored_rows .or. istep+nblk>na .or. (na/np_rows<=256 .and. nstor>=32)) then
+
+          ! Calculate scalar products of stored vectors.
+          ! This can be done in different ways, we use zherk
+
+          tmat = 0
+          if (l_rows>0) &
+             call PRECISION_HERK('U', 'C', nstor, l_rows, CONE, hvm, ubound(hvm,dim=1), CZERO, tmat, max_stored_rows)
+          nc = 0
+          do n=1,nstor-1
+            h1(nc+1:nc+n) = tmat(1:n,n+1)
+            nc = nc+n
+          enddo
+#ifdef WITH_MPI
+#ifdef HAVE_DETAILED_TIMINGS
+          call timer%start("mpi_communication")
+#endif
+          if (nc>0) call mpi_allreduce(h1, h2, nc, MPI_COMPLEX_PRECISION, MPI_SUM, mpi_comm_rows, mpierr)
+#ifdef HAVE_DETAILED_TIMINGS
+          call timer%stop("mpi_communication")
+#endif
+
+#else /* WITH_MPI */
+          if (nc>0) h2=h1
+#endif /* WITH_MPI */
+          ! Calculate triangular matrix T
+
+          nc = 0
+          tmat(1,1) = tau(ice-nstor+1)
+          do n=1,nstor-1
+            call PRECISION_TRMV('L', 'C', 'N', n, tmat, max_stored_rows, h2(nc+1),1)
+            tmat(n+1,1:n) = -conjg(h2(nc+1:nc+n))*tau(ice-nstor+n+1)
+            tmat(n+1,n+1) = tau(ice-nstor+n+1)
+            nc = nc+n
+          enddo
+
+          ! Q = Q - V * T * V**T * Q
+
+          if (l_rows>0) then
+            call PRECISION_GEMM('C', 'N', nstor, l_cols, l_rows, CONE, hvm, ubound(hvm,dim=1), &
+                        q, ldq, CZERO, tmp1 ,nstor)
+          else
+            tmp1(1:l_cols*nstor) = 0
+          endif
+#ifdef WITH_MPI
+#ifdef HAVE_DETAILED_TIMINGS
+          call timer%start("mpi_communication")
+#endif
+
+          call mpi_allreduce(tmp1, tmp2, nstor*l_cols, MPI_COMPLEX_PRECISION, MPI_SUM, mpi_comm_rows, mpierr)
+#ifdef HAVE_DETAILED_TIMINGS
+          call timer%stop("mpi_communication")
+#endif
+
+          if (l_rows>0) then
+            call PRECISION_TRMM('L', 'L', 'N', 'N', nstor, l_cols, CONE, tmat, max_stored_rows, tmp2, nstor)
+            call PRECISION_GEMM('N', 'N', l_rows, l_cols, nstor, -CONE, hvm, ubound(hvm,dim=1), &
+                        tmp2, nstor, CONE, q, ldq)
+          endif
+
+#else
+!          tmp2 = tmp1
+
+          if (l_rows>0) then
+            call PRECISION_TRMM('L', 'L', 'N', 'N', nstor, l_cols, CONE, tmat, max_stored_rows, tmp1, nstor)
+            call PRECISION_GEMM('N', 'N', l_rows, l_cols, nstor, -CONE, hvm, ubound(hvm,dim=1), &
+                        tmp1, nstor, CONE, q, ldq)
+          endif
+
+#endif
+
+!          if (l_rows>0) then
+!            call PRECISION_TRMM('L', 'L', 'N', 'N', nstor, l_cols, CONE, tmat, max_stored_rows, tmp2, nstor)
+!            call PRECISION_GEMM('N', 'N', l_rows, l_cols, nstor, -CONE, hvm, ubound(hvm,dim=1), &
+!                        tmp2, nstor, CONE, q, ldq)
+!          endif
+
+          nstor = 0
+        endif
+
+      enddo
+
+      deallocate(tmat, h1, h2, tmp1, tmp2, hvb, hvm, stat=istat, errmsg=errorMessage)
+      if (istat .ne. 0) then
+       print *,"trans_ev_complex: error when deallocating hvb "//errorMessage
+       stop
+      endif
+
+      call timer%stop("trans_ev_complex" // PRECISION_SUFFIX)
+
+    end subroutine trans_ev_complex_PRECISION