      program Oseen11

C---- Oseen11 calculates the 11 x 11 grand resistance matrix
c     of a rigid object build from spherical primary particles.
c     The primary spheres should not touch or overlap.
c
c     For details, consult
c       'Efficient Brownian Dynamics of rigid colloids
c         in linear flow fields based on the grand mobility matrix'
c       by
c       Duraivelan Palanisamy and Wouter K. den Otter
c       J. Chem. Phys. 148, 194112 (2018).
c
c     This program comes without guarantees.
c     Should you find an error, please let us know.
c     * Eq C8b is incorrect -- the correct expression is used here.
c
c     The 11x11 grand mobility matrix is not symmetric;
c     it will be after replacing -E by E in Eq. 10,
c     i.e. by multiplying the right 5 columns with -1.
c
c     By far the most time-consuming part of this code
c     is inverting the ( 11N x 11N ) matrix Mu.
c     In this release, both lapack and straightforward jacobi are used.
c     When substituting matInv with your favourite inversion routine,
c     note that all matrices to be inverted are symmetric.
c
c     The radius of the primary particle is the unit of length,
c     the viscosity is taken as 1.
c     Factors of 6 pi are included -- a single sphere yields Stokes.
c
c     input:
c     * in.config:
c       # title
c       # number of primary spheres
c       # absolute x, y and z coordinates, one sphere per line
c       
c     output:
c     * progress information to screen
c     * out.rstnc.abs
c         11 x 11 grand resistance matrix in absolute coordinates
c     * out.rstnc.rel
c         11 x 11 grand resistance matrix relative to hydro. center
c     * out.mblty.abs
c         11 x 11 grand mobility matrix in absolute coordinates
c     * out.mblty.rel
c         11 x 11 grand mobility matrix relative to hydro. center
c       The latter four are in ascii or binary (compile with -Dbin).
c
c     Compilation options:
c       -Dbin: binary output files
c       -Dlapack   use lapack / blas (don't forget to link libraries)
c                  i.s.o. jacobi for big matrix inversion
c       -mcmodel=medium for large arrays
c     e.g.
c       ifort -cpp -O2 -xhost -Dlapack -Dbin -mkl \\
c             -mcmodel=medium -shared-intel \\
c             -o Oseen11 Oseen11.f
c       gfortran -O2 -cpp -Dlapack -Dbin -llapack -mcmodel=medium \\
c             -o Oseen11 Oseen11.f
c
c     Please cite the above manuscript if you use this code.
c
c     For comments and suggestions: w.k.denotter@utwente.nl
c
c     Enschede, 04 June 2020


      implicit none

      integer mxN, mxN11
      parameter ( mxN   =    2 500 )
      parameter ( mxN11 = mxN * 11 )

      real*8  r( 3, mxN )
      integer N

      integer i, j, k, l
      integer ii, jj, kk, ll
      integer N11
      integer a, b, c, d
      integer aa, bb, cc, dd
      real*8  rija( 3 ), rijab( 3, 3 ), dijab( 3, 3 ), pijab( 3, 3 )
      real*8  Kron( 3, 3 )
      real*8  Levi( 3, 3, 3 )
      real*8  Mvf( 3, 3 ), Mwf( 3, 3 ), Mwt( 3, 3 )
      real*8  Mef( 3, 3, 3 ), Met( 3, 3, 3 ), Mes( 3, 3, 3, 3 )

      real*8  Mu ( mxN11, mxN11 )
      real*8  Xi ( mxN11, mxN11 )
      real*8  rst( 11, 11 ), mbl( 11, 11 )
      real*8  tmp( 11, 11 ), mblOld( 11, 11 )

      real*8  Sv2m( 3, 3, 5 )
      real*8  Sm2v( 5, 3, 3 ) ! note: order flipped for readability code
      real*8  Ev2m( 3, 3, 5 )
      real*8  Em2v( 5, 3, 3 ) ! note: order flipped for readability code

      real*8  eps( 3, 3, mxN )
      real*8  phi( 3, 5, mxN )
      real*8  psi( 5, 3, mxN )
      real*8  epsDx( 3, 3 ), phiDx( 3, 5 ), psiDx( 5, 3 )

      real*8  rij, rij1, rij2, rij3, rij4, rij5
      real*8  dij, jid
      real*8  xvf, yvf, ywf, xwt, ywt, xef, yef, yet, xes, yes, zes
      real*8  minr
      real*8  sm, s, mx

      real*8  cntr( 3 )
      real*8  mat( 3, 3 ), tam( 3, 3 ), vec( 3 )

      integer iia, iiaa, jjb, jjbb, jjc, jjcc, iid
      integer iik, jjl

      character*80 title

      real*8  pi

#ifdef lapack

      integer lwork
      parameter ( lwork = mxN11 * 64 )

      integer ipiv( mxN11 ), info
      real*8  work( lwork )

      character uplo
      parameter ( uplo = 'U' )

      real*8  evc( 3, 3 ), evl( 3 )

#else

      real*8  Mu2( mxN11, mxN11 )
      real*8  evc( mxN11, mxN11 ), evl( mxN11 )

#endif

C----

      print *
      print *,'   ###############'
      print *,'   #             #'
      print *,'   #   Oseen11   #'
      print *,'   #             #'
      print *,'   ###############'
      print *
      print *,'Authors: Duraivelan Palanisamy and Wouter K. den Otter.'
      print *,'All rights reserved.'
      print *,'Enschede, 04 June 2020.'
      print *
      print *,'This code comes with no warrenties whatsoever.'
      print *
      print *,'For the underlying theory, consult the following paper:'
      print *,'  Efficient Brownian Dynamics of rigid colloids',
     x           ' in linear flow fields'
      print *,'      based on the grand mobility matrix,'
      print *,'    Duraivelan Palanisamy and Wouter K. den Otter,'
      print *,'    J. Chem. Phys. 148, 194112 (2018).'
      print *
      print *,'If you use this code, please cite above paper.'
      print *
      print *,'For comments and suggestions: w.k.denotter@utwente.nl'
      print *

      open( 20, file = 'in.config',     form   = 'formatted',
     x                                  status = 'old' )

#ifdef bin
      open( 30, file = 'out.rstnc.abs', form   = 'unformatted',
     x                                  status = 'new' )
      open( 31, file = 'out.mblty.abs', form   = 'unformatted',
     x                                  status = 'new' )
      open( 32, file = 'out.mblty.rel', form   = 'unformatted',
     x                                  status = 'new' )
      open( 33, file = 'out.rstnc.rel', form   = 'unformatted',
     x                                  status = 'new' )
#else
      open( 30, file = 'out.rstnc.abs', form   = 'formatted',
     x                                  status = 'new' )
      open( 31, file = 'out.mblty.abs', form   = 'formatted',
     x                                  status = 'new' )
      open( 32, file = 'out.mblty.rel', form   = 'formatted',
     x                                  status = 'new' )
      open( 33, file = 'out.rstnc.rel', form   = 'formatted',
     x                                  status = 'new' )
#endif

      print *, '---- reading configuration file'
      print *

      read( 20, '(a80)' ) title
      read( 20, * ) N
      print *, title
      print *, N
      print *
      if ( N .gt. mxN )
     x  stop 'N > mxN; increase array size and recompile'
      do i = 1, N
        read( 20, * ) ( r( j, i ), j = 1, 3 )
      enddo ! i
      close( 20 )

!     N = 10
!     print *,'***** reduced to ten spheres for testing ******'
!     print *

#ifdef bin
      write( 30 ) title
      write( 31 ) title
      write( 32 ) title
      write( 33 ) title
#else
      write( 30, * ) title
      write( 31, * ) title
      write( 32, * ) title
      write( 33, * ) title
#endif

      N11 = 11 * N

      Kron(  :,1) = [  1,  0,  0 ]
      Kron(  :,2) = [  0,  1,  0 ]
      Kron(  :,3) = [  0,  0,  1 ]

      Levi(1,:,1) = [  0,  0,  0 ]
      Levi(2,:,1) = [  0,  0,  1 ]
      Levi(3,:,1) = [  0, -1,  0 ]
      Levi(1,:,2) = [  0,  0, -1 ]
      Levi(2,:,2) = [  0,  0,  0 ]
      Levi(3,:,2) = [  1,  0,  0 ]
      Levi(1,:,3) = [  0,  1,  0 ]
      Levi(2,:,3) = [ -1,  0,  0 ]
      Levi(3,:,3) = [  0,  0,  0 ]

C---- vector to matrix and reverse, for stress and strain

      Sv2m(1,:,1) = [ 0, 1, 0 ]  !  Eq B2
      Sv2m(2,:,1) = [ 1, 0, 0 ]
      Sv2m(3,:,1) = [ 0, 0, 0 ]

      Sv2m(1,:,2) = [ 0, 0, 1 ]
      Sv2m(2,:,2) = [ 0, 0, 0 ]
      Sv2m(3,:,2) = [ 1, 0, 0 ]

      Sv2m(1,:,3) = [ 0, 0, 0 ]
      Sv2m(2,:,3) = [ 0, 0, 1 ]
      Sv2m(3,:,3) = [ 0, 1, 0 ]

      Sv2m(1,:,4) = [ 2,  0,  0 ] / 3.
      Sv2m(2,:,4) = [ 0, -1,  0 ] / 3.
      Sv2m(3,:,4) = [ 0,  0, -1 ] / 3.

      Sv2m(1,:,5) = [ 1,  0,   0 ] / 3.
      Sv2m(2,:,5) = [ 0,  1,   0 ] / 3.
      Sv2m(3,:,5) = [ 0,  0,  -2 ] / 3.

      Sm2v(1,:,:) = Sv2m(:,:,1) / 2.  !  Eq B3
      Sm2v(2,:,:) = Sv2m(:,:,2) / 2.
      Sm2v(3,:,:) = Sv2m(:,:,3) / 2.

      Sm2v(4,1,:) = [ 1,  0,  0 ]
      Sm2v(4,2,:) = [ 0, -1,  0 ]
      Sm2v(4,3,:) = [ 0,  0,  0 ]

      Sm2v(5,1,:) = [ 0,  0,   0 ]
      Sm2v(5,2,:) = [ 0,  1,   0 ]
      Sm2v(5,3,:) = [ 0,  0,  -1 ]

      do i = 1, 5
        Ev2m( :, :, i ) = Sm2v( i, :, : )  !  Eq B4
        Em2v( i, :, : ) = Sv2m( :, :, i )
      enddo ! i

      pi = 4. * atan( 1. )

      minr  = 9999.

C---- compile matrix ( u, w, -e ) <-- ( f, t, s )
c     by generating the N x N ( 11 x 11 ) i-j submatrices.

      print *, '---- filling ( 11N x 11N ) mobility matrix'
      print *

      Mu = 0.

      do i = 1, N
        
        do j = 1, N  !  can be improved, but not worth the effort

C----     compose i-j matrices

          if ( i .eq. j ) then
            dij  = 1  !  Kronecker
            jid  = 0  !  1 - Kronecker
            rija = 0.
            rij  = 0.
            rij1 = 0.
            rij2 = 0.
            rij3 = 0.
            rij4 = 0.
            rij5 = 0.
          else
            dij  = 0
            jid  = 1
            rija = r( :, j ) - r( :, i )
            rij  = sqrt( rija( 1 ) **2 + rija( 2 ) **2 + rija( 3 ) **2 )
            minr = min( minr, rij )
            rij1 = 1. / rij
            rij2 =      rij1 ** 2
            rij3 =      rij1 ** 3
            rij4 =      rij1 ** 4
            rij5 =      rij1 ** 5
            rija = rija * rij1
          endif

C----     distance dependence, Eq A5

          xvf =        dij + jid * ( 1.5   * rij1 -       rij3 )
          yvf =        dij + jid * ( 0.75  * rij1 + 0.5 * rij3 )
          ywf =            - jid *   0.75  * rij2
          xwt = 0.75 * dij + jid *   0.75  * rij3
          ywt = 0.75 * dij - jid *   0.375 * rij3
          xef =              jid * ( 2.25  * rij2 - 3.6 * rij4 )
          yef =              jid *   1.2   * rij4
          yet =            - jid *   1.125 * rij3
          xes = 1.35 * dij - jid * ( 6.75  * rij3 - 16.2 * rij5 )
          yes = 0.45 * dij + jid * ( 1.125 * rij3 -  3.6 * rij5 )
          zes = 0.45 * dij + jid *   0.9   * rij5

C----     basis tensors

          do a = 1, 3
            do b = 1, 3
              rijab( a, b ) = rija( a ) * rija( b )
            enddo ! j
          enddo ! i
          dijab = rijab - Kron / 3.
          pijab = Kron - rijab

C----     combination of the above, Eq A3

          Mvf = xvf * rijab + yvf * pijab

          Mwf = 0.
          do a = 1, 3
            do b = 1, 3
              do c = 1, 3
               Mwf( a, b ) = Mwf( a, b )
     x              + ywf * Levi( a, b, c ) * rija( c )
              enddo ! c
            enddo ! b
          enddo ! a

          Mwt = xwt * rijab + ywt * pijab

          Mef = 0.
          do a = 1, 3
            do b = 1, 3
              do c = 1, 3
                Mef( a, b, c ) = Mef( a, b, c )
     x              + xef * dijab( a, b ) * rija( c )
     x              + yef * ( rija( a ) * pijab( b, c )
     x                      + rija( b ) * pijab( a, c ) )
              enddo ! c
            enddo ! b
          enddo ! a

          Met = 0.
          do a = 1, 3
            do b = 1, 3
              do c = 1, 3
                do d = 1, 3
                  Met( a, b, c ) = Met( a, b, c )
     x                + yet * ( rija( a ) * Levi( b, c, d )
     x                        + rija( b ) * Levi( a, c, d ) )
     x                                          * rija( d )
                enddo ! d
              enddo ! c
            enddo ! b
          enddo ! a

          Mes = 0.
          do a = 1, 3
            do b = 1, 3
              do c = 1, 3
                do d = 1, 3
                  Mes( a, b, c, d ) = Mes( a, b, c, d )
     x              + xes *   dijab( a, b ) * dijab( c, d )
     x              + yes * ( rija( a ) * pijab( b, c ) * rija( d )
     x                      + rija( a ) * pijab( b, d ) * rija( c )
     x                      + rija( b ) * pijab( a, c ) * rija( d )
     x                      + rija( b ) * pijab( a, d ) * rija( c ) )
     x              + zes * ( pijab( a, c ) * pijab( b, d )
     x                      + pijab( b, c ) * pijab( a, d )
     x                      - pijab( a, b ) * pijab( c, d ) )
                enddo ! d
              enddo ! c
            enddo ! b
          enddo ! a

C----     insert i-j matrices in 11N x 11N matrix Mu

          ii = 11 * ( i - 1 )
          jj = 11 * ( j - 1 )

          do a = 1, 3
            iia  = ii  + a
            iiaa = iia + 3
            do b = 1, 3
              jjb  = jj  + b
              jjbb = jjb + 3
              Mu( iia , jjb  ) = Mvf( a, b )
              Mu( iiaa, jjb  ) = Mwf( a, b )
              Mu( iiaa, jjbb ) = Mwt( a, b )
              Mu( jjb , iiaa ) = Mwf( a, b )  !  Eq A4a
            enddo ! b
          enddo ! a

          do c = 1, 3
            jjc  = jj  + c
            jjcc = jjc + 3
            do a = 1, 3
              do b = 1, 3
                do k = 1, 5  !  k <-- ( a, b )
                  iik = ii + k + 6
                  Mu( iik, jjc  ) = Mu  ( iik,       jjc  )
     x                            + Em2v(   k, a, b       )
     x                            * Mef (      a, b,   c  )
                  Mu( iik, jjcc ) = Mu  ( iik,       jjcc )
     x                            + Em2v(   k, a, b       )
     x                            * Met (      a, b,   c  )
                enddo ! k
              enddo ! b
            enddo ! a
            do k = 1, 5
              iik = ii + k + 6
              Mu( jjc,  iik ) = Mu( iik, jjc  )  !  Eq A4b
              Mu( jjcc, iik ) = Mu( iik, jjcc )  !  Eq A4c
            enddo ! k
          enddo ! c

          do a = 1, 3
            do b = 1, 3
              do k = 1, 5  !  k <-- ( a, b )
                iik = ii + k + 6
                do c = 1, 3
                  do d = 1, 3
                    do l = 1, 5  !  l <-- ( c, d )
                      jjl = jj + l + 6
                      Mu( iik, jjl ) =  Mu  ( iik,             jjl )
     x                                + Em2v(   k, a, b            )
     x                                * Mes (      a, b, c, d      )
     x                                * Sv2m(            c, d,   l )
                    enddo ! l
                  enddo ! d
                enddo ! c
              enddo ! k
            enddo ! b
          enddo ! a

        enddo ! j
      enddo ! i

      print '(a20,f10.5)', 'minimum distance:', minr
      print *
      if ( minr .le. 2. ) stop 'minimum distance too low'

C---- symmetry test

      sm = 0.
      mx = 0.
      do i = 1, N11
        do j = i + 1, N11
          s  = abs( Mu( i, j ) - Mu( j, i ) )
          sm = sm + s
          mx = max( mx, s )
        enddo ! j
      enddo ! i
      print '(a20,2e12.4)', 'symmetry test:',
     x       sm / dble( N11 ** 2 ), mx
      print *

      Mu  = Mu / ( 6. * pi )

C---- invert the ( 11N x 11N ) grand mobility matrix

      print *, '---- inverting mobility matrix (be patient...)'
      print *

#ifdef lapack

      Xi = Mu

c     ignoring symmetry

!     call dgetrf( N11, N11, Xi, mxN11, ipiv, info )
!     if ( info .ne. 0 ) stop 'mobility matrix is singular'
!     call dgetri( N11, Xi, mxN11, ipiv, work, lwork, info )
!     if ( info .ne. 0 ) stop 'inversion failed'

c     taking advantage of symmetry

      if ( uplo .ne. 'U' ) stop 'This code assumes uplo = U'

      call dsytrf( uplo, N11, Xi, mxN11, ipiv, work, lwork, info)

      if ( info .ne. 0 ) stop 'problem in dsytrf'

      call dsytri( uplo, N11, Xi, mxN11, ipiv, work, info )

      if ( info .ne. 0 ) stop 'inversion failed'

      do i = 1, N11 - 1  !  this is based on uplo = U.
        do j = i + 1, N11
          Xi( j, i ) = Xi( i, j )
        enddo ! j
      enddo ! i

#else

      Mu2 = Mu

      call matInv( Mu2, N11, mxN11, Xi, evl, evc )  !  Mu2 is destroyed

#endif

      print *,'Section of inverse'
      do i = 1, 11
        print '(11f15.5)', ( Xi( i, j ), j = 1, 11 )
      enddo ! i
      print *

      sm = 0.
      do i = 1, N11
        do j = 1, N11
          if ( i .eq. j ) then
            s = - 1.
          else
            s = 0.
          endif
          do k = 1, N11
            s = s + Mu( i, k ) * Xi( k, j )
          enddo ! k
        sm = sm + s ** 2
        enddo ! j
      enddo ! i
      print '(a20,e12.4)', 'inversion test:', sm
      print *

C---- preparation of auxilliary matrices eps, phi and psi

      print *, '---- reduction to ( 11 x 11 ) resistance matrix'
      print *

      eps = 0.
      phi = 0.
      psi = 0.

      do i = 1, N

        do a = 1, 3
          do b = 1, 3
            do c = 1, 3
              eps( a, b, i ) = eps ( a,    b, i )
     x                       + Levi( a, c, b    )
     x                       *    r(    c,    i )
            enddo ! c
          enddo ! b
        enddo ! b

        do a = 1, 3
          do k = 1, 5
            do b = 1, 3
              phi( a, k, i ) = phi ( a,    k, i )
     x                       + Ev2m( a, b, k    )
     x                       *    r(    b,    i )
              psi( k, a, i ) = psi ( k,    a, i )
     x                       + Sm2v( k, b, a    )
     x                       *    r(    b,    i )
            enddo ! b
          enddo ! k
        enddo ! a

      enddo ! i

C---- reduction to a ( 11 x 11 ) grand mobility matrix

      rst = 0.

      do i = 1, N
        ii = 11 * ( i - 1 )

        do j = 1, N
          jj = 11 * ( j - 1 )

C----     the ( 3 x 3 ) sub matrices: Eqs C4a, C4b, C5a and C5b

          do a = 1, 3
              aa =  a + 3
            iia  = ii +  a    !  f,i,a
            iiaa = ii + aa    !  t,i,a

            do b = 1, 3
                bb =  b + 3
              jjb  = jj +  b  !  v,j,b
              jjbb = jj + bb  !  w,j,b

              rst(  a,  b ) = rst(  a,  b ) + Xi( iia,  jjb  ) !  Xfv
              rst(  a, bb ) = rst(  a, bb ) + Xi( iia,  jjbb ) !  Xfw
              rst( aa,  b ) = rst( aa,  b ) + Xi( iiaa, jjb  ) !  Xtv
              rst( aa, bb ) = rst( aa, bb ) + Xi( iiaa, jjbb ) !  Xtw

              do c = 1, 3  !  Xfw - Xfv eps
                jjc = jj + c
                rst( a, bb ) = rst(   a,                 bb    )
     x                       - Xi ( iia, jjc ) * eps( c,  b, j )
              enddo ! c

              do d = 1, 3  !  Xtv + eps Xfv
                iid = ii + d
                rst( aa, b ) = rst( aa,                     b )
     x                       + eps(  a, d, i ) * Xi( iid, jjb )
              enddo ! d

              do c = 1, 3  !  Xtw - Xtv eps
                jjc = jj + c
                rst( aa, bb ) = rst(   aa,                 bb     )
     x                        - Xi ( iiaa, jjc ) * eps( c,  b , j )
              enddo ! c
              do d = 1, 3  !      + eps Xfw
                iid  = ii + d
                rst( aa, bb ) = rst( aa,                      bb )
     x                        + eps(  a, d, i ) * Xi ( iid, jjbb )
              enddo ! d
              do c = 1, 3  !      - eps Xfv eps
                jjc = jj + c
                do d = 1, 3
                  iid  = ii + d
                  rst( aa, bb ) = rst( aa,           bb    )
     x                          - eps(  a,   d,          i )
     x                          * Xi (     iid, jjc        )
     x                          * eps(            c,  b, j )
                enddo ! d
              enddo ! c

            enddo ! b

          enddo ! a

C----     the ( 3 x 5 ) sub matrices, Eqs C4c and C5c

          do a = 1, 3
              aa =  a +  3
            iia  = ii +  a    !  f,i,a
            iiaa = ii + aa    !  t,i,a

            do l = 1, 5
                ll =  l + 6
              jjl  = jj + ll  !  e,j,l

              rst(  a, ll ) = rst(  a, ll ) + Xi( iia,  jjl )  !  Xfe
              rst( aa, ll ) = rst( aa, ll ) + Xi( iiaa, jjl )  !  Xte

              do c = 1, 3  !  Xfe + Xfv phi
                jjc = jj + c
                rst( a, ll ) = rst(   a,                 ll    )
     x                       + Xi ( iia, jjc ) * phi( c,  l, j )
              enddo ! c

              do c = 1, 3  !  Xte + Xtv phi
                jjc = jj + c
                rst( aa, ll ) = rst(   aa,                 ll    )
     x                        + Xi ( iiaa, jjc ) * phi( c,  l, j )
              enddo ! c
              do d = 1, 3  !      + eps Xfe
                iid  = ii + d
                rst( aa, ll ) = rst( aa,                      ll )
     x                        + eps(  a, d, i ) * Xi ( iid, jjl  )
              enddo ! d
              do c = 1, 3  !      + eps Xfv phi
                jjc = jj + c
                do d = 1, 3
                  iid  = ii + d
                  rst( aa, ll ) = rst( aa,           ll    )
     x                          + eps(  a,   d,          i )
     x                          * Xi (     iid, jjc        )
     x                          * phi(            c,  l, j )
                enddo ! d
              enddo ! c

            enddo ! l

          enddo ! a
          
C----     the ( 5 x 3 ) sub matrices, Eqs C8a and C8b

          do k = 1, 5
              kk =  k + 6
            iik  = ii + kk    !  s,i,k

            do b = 1, 3
                bb =  b + 3
              jjb  = jj +  b  !  v,j,b
              jjbb = jj + bb  !  w,j,b

              rst( kk,  b ) = rst( kk,  b ) + Xi( iik, jjb  )  !  Xsv
              rst( kk, bb ) = rst( kk, bb ) + Xi( iik, jjbb )  !  Xsw

              do d = 1, 3  !  Xsv + psi Xfv
                iid = ii + d
                rst( kk, b ) = rst( kk,                     b )
     x                       + psi(  k, d, i ) * Xi( iid, jjb )
              enddo ! d

              do c = 1, 3  !  Xsw - Xsv eps      !  wrong in eq C8b
                jjc = jj + c
                rst( kk, bb ) = rst(   kk,                 bb    )
     x                        - Xi ( iik,  jjc ) * eps( c,  b, j )
              enddo ! c
              do d = 1, 3  !      + psi Xfw
                iid = ii + d
                rst( kk, bb ) = rst( kk,                     bb )
     x                        + psi(  k, d, i ) * Xi( iid, jjbb )
              enddo ! d
              do c = 1, 3  !      - psi Xfv eps  !  wrong in Eq C8b
                jjc = jj + c
                do d = 1, 3
                  iid = ii + d
                  rst( kk, bb ) = rst( kk,           bb    )
     x                          - psi(  k,   d,          i )
     x                          * Xi (     iid, jjc        )
     x                          * eps(            c,  b, j )
                enddo ! d
              enddo ! c

            enddo ! a

          enddo ! k

C----     the ( 5 x 5 ) sub matrix, Eq. C8c

          do k = 1, 5
              kk =  k + 6
            iik  = ii + kk    !  s,i,k

            do l = 1, 5
                ll =  l + 6
              jjl  = jj + ll  !  e,j,l

              rst( kk, ll ) = rst( kk, ll ) + Xi( iik, jjl )  !  Xse

              do c = 1, 3  !  Xse + Xsv phi
                jjc = jj + c
                rst( kk, ll ) = rst(   kk,                 ll    )
     x                        + Xi ( iik,  jjc ) * phi( c,  l, j )
              enddo ! c
              do d = 1, 3  !      + psi Xfe
                iid = ii + d
                rst( kk, ll ) = rst( kk,                     ll )
     x                        + psi(  k, d, i ) * Xi( iid, jjl  )
              enddo ! d
              do c = 1, 3  !      + psi Xfv phi
                jjc = jj + c
                do d = 1, 3
                  iid = ii + d
                  rst( kk, ll ) = rst( kk,           ll    )
     x                          + psi(  k,   d,          i )
     x                          * Xi (     iid, jjc        )
     x                          * phi(            c,  l, j )
                enddo ! d
              enddo ! c

            enddo ! l

          enddo ! k

        enddo ! j
      enddo ! i

C---- ouput

      print *,'---- in absolute coordinates'
      print *

      print *,'grand resistance matrix'
      print '(11f15.5)', rst
      print *

#ifdef bin
      write( 30 ) rst
#else
      do i = 1, 11
        write( 30, * ) i
        do j = 1, 11
          write( 30, * ) rst( j, i )
        enddo ! j
      enddo ! i
#endif
      close( 30 )

C---- symmetry test

      sm = 0.
      mx = 0.
      do i = 1, 11
        do j = i + 1, 11
          s  = abs( rst( i, j ) - rst( j, i ) )
          sm = sm + s
          mx = max( mx, s )
        enddo ! j
      enddo ! i
      print '(a20,2e12.4)', 'symmetry test:', sm / 121., mx
      print *

C---- determine the mobility matrix

      tmp = rst

      call partInv( rst, mbl )

      print *,'grand mobility matrix'
      print '(11f15.5)', mbl
      print *

#ifdef bin
      write( 31 ) mbl
#else
      do i = 1, 11
        write( 31, * ) i
        do j = 1, 11
          write( 31, * ) mbl( j, i )
        enddo ! j
      enddo ! i
#endif
      close( 31 )

C---- symmetry test

      sm = 0.
      mx = 0.
      mblOld = mbl
      mblOld( :, 7 : 11 ) = - mblOld( :, 7 : 11 )
      do i = 1, 11
        do j = i + 1, 11
          s  = abs( mblOld( i, j ) - mblOld( j, i ) )
          sm = sm + s
          mx = max( mx, s )
        enddo ! j
      enddo ! i
      print '(a20,2e12.4)', 'pseudo-symmetry:', sm / 121., mx
      print *

C---- deviation from hydrodynamic center

      sm = 0.
      do a = 1, 3
        do b = 1, 3
          sm = sm + ( mbl( a + 3, b ) - mbl( b + 3, a ) ) ** 2
        enddo ! b
      enddo ! a
      print '(a20,e12.4)', 'off-center:', sm
      print *

C---- hydrodynamic center: mobility based and in current coordinates

      cntr = 0  !  Eq G2
      vec  = 0.
      do a = 1, 3
        do b = 1, 3
          do c = 1, 3
            cc = c + 3
            vec( a ) = vec( a ) + Levi( a, b,  c )
     x                          * mbl (    b, cc )
          enddo ! c
        enddo ! b
      enddo ! a
      mat = mbl( 4 : 6, 4 : 6 )
      sm  = mat( 1, 1 ) + mat( 2, 2 ) + mat( 3, 3 )
      mat = mat - sm * Kron

      call matInv( mat, 3, 3, tam, evl, evc )  !  Eq. G2

      do a = 1, 3
        do b = 1, 3
          cntr( a ) = cntr( a ) + tam( a, b ) * vec( b )
        enddo ! b
      enddo ! a

      print *,'location mobility center'
      print '(3f15.5)', cntr
      print *

C---- mobility matrix relative to mobility center:
C     subtract cntr from all coordinates and repeat above lines
c     or apply a simple transformation.

      epsDx = 0.
      do a = 1, 3
        do b = 1, 3
          do c = 1, 3
            epsDx( a, b ) = epsDx( a,    b )
     x                    + Levi ( a, c, b )
     x                    *  cntr(    c    )
          enddo ! c
        enddo ! b
      enddo ! b

      phiDx = 0.
      psiDx = 0.
      do a = 1, 3
        do k = 1, 5
          do b = 1, 3
            phiDx( a, k ) = phiDx( a,    k )
     x                    + Ev2m ( a, b, k )
     x                    *  cntr(    b    )
            psiDx( k, a ) = psiDx( k,    a )
     x                    + Sm2v ( k, b, a )
     x                    *  cntr(    b    )
          enddo ! b
        enddo ! k
      enddo ! a

      mblOld = mbl

C---- the ( 3 x 3 ) sub matrices: Eqs G1a, G1b, G1d

      do a = 1, 3      !  v
        aa = a + 3     !  w
        do b = 1, 3    !  f
          bb = b + 3   !  t
          do c = 1, 3
            cc = c + 3
            mbl(  a,  b ) = mbl   (  a,                     b )
     x                    + mblOld(  a, cc ) * epsDx (  c,  b )
     x                    - epsDx (  a,  c ) * mblOld( cc,  b )
            mbl(  a, bb ) = mbl   (  a,                    bb )
     x                    - epsDx (  a, c )  * mblOld( cc, bb )
            mbl( aa,  b ) = mbl   ( aa,                     b )
     x                    + mblOld( aa, cc ) * epsDx (  c,  b )
            do d = 1, 3
              dd = d + 3
              mbl( a, b ) = mbl   ( a,         b )
     x                    - epsDx ( a,  c        )
     x                    * mblOld(    cc, dd    )
     x                    * epsDx (         d, b )
            enddo ! d
          enddo ! c
        enddo ! b
      enddo ! a

C---- the ( 3 x 5 ) and ( 5 x 3 ) sub matrices: Eq G1c and G1e

      do a = 1, 3     !  v, f
        do k = 1, 5
          kk = k + 6  !  e, s
          do c = 1, 3
            cc = c + 3
            mbl(  a, kk ) = mbl   (  a,                    kk ) 
     x                    - epsDx (  a,  c ) * mblOld( cc, kk )
            mbl( kk,  a ) = mbl   ( kk,                     a ) 
     x                    + mblOld( kk, cc ) * epsDx (  c,  a )
          enddo ! c
        enddo ! k
      enddo ! a
      mbl( 1 :  3, 7 : 11 ) = mbl( 1 :  3, 7 : 11 ) + phiDx
      mbl( 7 : 11, 1 :  3 ) = mbl( 7 : 11, 1 :  3 ) - psiDx

C---- output

      print *,'---- in relative coordinates'
      print *

      print *,'grand mobility matrix'
      print '(11f15.5)', mbl
      print *

#ifdef bin
      write( 32 ) mbl
#else
      do i = 1, 11
        write( 32, * ) i
        do j = 1, 11
          write( 32, * ) mbl( j, i )
        enddo ! j
      enddo ! i
#endif
      close( 32 )

C---- symmetry test

      sm = 0.
      mx = 0.
      mblOld = mbl
      mblOld( :, 7 : 11 ) = - mblOld( :, 7 : 11 )
      do i = 1, 11
        do j = i + 1, 11
          s  = abs( mblOld( i, j ) - mblOld( j, i ) )
          sm = sm + s
          mx = max( mx, s )
        enddo ! j
      enddo ! i
      print '(a20,2e12.4)', 'pseudo-symmetry:', sm / 121., mx
      print *

C---- deviation from hydrodynamic center

      sm = 0.
      do a = 1, 3
        do b = 1, 3
          sm = sm + ( mbl( a + 3, b ) - mbl( b + 3, a ) ) ** 2
        enddo ! b
      enddo ! a
      print '(a20,e12.4)', 'off-center:', sm
      print *

C---- determine the resistance matrix

      tmp = mbl

      call partInv( mbl, rst )

      print *,'grand resistance matrix'
      print '(11f15.5)', rst
      print *

#ifdef bin
      write( 33 ) rst
#else
      do i = 1, 11
        write( 33, * ) i
        do j = 1, 11
          write( 33, * ) rst( j, i )
        enddo ! j
      enddo ! i
#endif
      close( 33 )

C---- symmetry test

      sm = 0.
      mx = 0.
      do i = 1, 11
        do j = i + 1, 11
          s  = abs( rst( i, j ) - rst( j, i ) )
          sm = sm + s
          mx = max( mx, s )
        enddo ! j
      enddo ! i
      print '(a20,2e12.4)', 'symmetry test:', sm / 121., mx
      print *

      print *,'---- bye'
      print *

      end ! Oseen11

C======================================================================

      subroutine matInv( a, nn, np, b, evl, evc )

c     input:
c       a   = real symmetric matrix, upper triangle is destroyed
c       nn  = size of matrices
c       np  = size of arrays
c     output:
c       b    = inverse matrix
c       evl  = eigenvalues
c       evc  = normalized eigenvectors as columns

      implicit none

      integer nn, np
      real*8  a( np, np ), b( np, np ), evl( np ), evc( np, np )

      real*8  small
      parameter ( small = 1.e-10 )

      integer i, j, k

      call jacobi( a, nn, np, evl, evc, i )  !  a is destroyed...

      do i = 1, nn
        if ( abs( evl( i ) ) .le. small ) then
          print *, 'small eigenvalue'
          print *, i, evl( i )
          stop 'matInv stopped'
        endif ! evl
        evl( i ) = 1. / evl( i )
      enddo ! i

!     print *, ( evl( i ), i = 1, nn )

      b = 0.  
      do k = 1, nn
        do j = 1, nn
          do i = 1, nn
            b( i, j ) = b  ( i,                       j    )   
     x                + evc( i, k ) * evl( k ) * evc( j, k ) 
          enddo ! i 
        enddo ! j 
      enddo ! k 

      end ! matInv

C======================================================================

      subroutine partInv( Rin, Mout )

      implicit none

      real*8 Rin( 11, 11 ), Mout( 11, 11 )

      real*8 Q( 6, 6 ), Qinv( 6, 6 )
      real*8 R( 6, 5 ), S( 5, 6 ), T( 5, 5 )

      real*8 evc( 6, 6 ), evl( 6 )

      integer i, j, k, l
      integer ii, jj

C----

      Mout = 0.

      Q( 1 : 6, 1 : 6 ) = Rin( 1 :  6,  1 :  6 )
      R( 1 : 6, 1 : 5 ) = Rin( 1 :  6,  7 : 11 )
      S( 1 : 5, 1 : 6 ) = Rin( 7 : 11,  1 :  6 )
      T( 1 : 5, 1 : 5 ) = Rin( 7 : 11,  7 : 11 )

C---- top left, Q^-1

      call matInv( Q, 6, 6, Qinv, evl, evc )  !  Q is destroyed...

      Mout( 1:6, 1:6 ) = Qinv

C---- top right, - Q^-1 * R

      do i = 1, 6
        do j = 1, 5
          jj = j + 6
          do k = 1, 6
            Mout( i, jj ) = Mout( i,             jj )
     x                    - Qinv( i, k ) * R( k,  j )
          enddo ! k
        enddo ! j
      enddo ! i

C---- bottom left, S * Q^-1

      do i = 1, 5
        ii = i + 6
        do j = 1, 6
          do k = 1, 6
            Mout( ii, j ) = Mout( ii,                j )
     x                      + S (  i, k ) * Qinv( k, j )
          enddo ! k
        enddo ! j
      enddo ! i

C---- bottom right, T - S * Q^-1 R

      Mout( 7 : 11, 7 : 11 ) = T

      do i = 1, 5
        ii = i + 6
        do j = 1, 5
          jj = j + 6
          do k = 1, 6
            do l = 1, 6
              Mout( ii, jj )
     x             = Mout( ii,                            jj )
     x               - S (  i, k ) * Qinv( k, l ) * R( l,  j )
            enddo ! l
          enddo ! k
        enddo ! j
      enddo ! i

      end ! partInv

C======================================================================

      subroutine jacobi( a, nn, np, d, v, m )

c     Jacobi's approach to diagonalising a real symmetric matrix
c
c     Source:
c     numerical recipes in fortran 77,
c     W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
c     version 2.08, Cambridge University Press, New York, NY, U.S.A.,
c     section 11.1
c
c     input:
c       a  = real symmetric matrix, upper triangle is destroyed
c       nn = size of matrices
c       np = size of arrays
c     output:
c       d  = eigenvalues
c       v  = normalized eigenvectors as columns
c       m  = number of iterations

      implicit none

      integer nn, np, m
      real*8  a( np,  np ), d( np ), v( np, np )
      
#ifdef lapack
      integer     nmax
      parameter ( nmax =  100 )
#else
      integer     nmax
      parameter ( nmax =  27 500 )
#endif

      integer i, ip, iq, j
      real*8  c, g, h, s, sm, t, tau, theta, tresh
      real*8  b( nmax ),z( nmax )

      if ( nn .gt. nmax ) stop 'Jacobi: nn > nmax'

      do ip = 1, nn
        do iq = 1, nn
          v( iq, ip ) = 0.
        enddo ! iq
        v( ip, ip ) = 1.
        b( ip ) = a( ip, ip )
        d( ip ) = b( ip )
        z( ip ) = 0.
      enddo ! ip

      m = 0
      do i = 1, 50
        sm = 0.
        do ip = 1, nn - 1
          do iq = ip + 1, nn
            sm = sm + abs( a( ip, iq ) )
          enddo ! iq
        enddo ! ip
        if ( sm .eq. 0. ) return
        if ( i .lt. 4 ) then
          tresh = 0.2 * sm / ( nn ** 2 )
        else
          tresh = 0.
        endif ! i
        do ip = 1, nn - 1
          do iq = ip + 1, nn
            g = 100. * abs( a( ip, iq ) )
            if ( ( i .gt. 4 ) .and.
     x           ( abs( d( ip ) ) + g .eq. abs( d( ip ) ) ) .and.
     x           ( abs( d( iq ) ) + g .eq. abs( d( iq ) ) ) ) then
              a( ip, iq ) = 0.
            else if ( abs( a( ip, iq ) ) .gt. tresh ) then
              h = d( iq ) - d( ip )
              if ( abs( h ) + g .eq. abs( h ) ) then
                t = a( ip, iq ) / h
              else
                theta = 0.5 * h / a( ip, iq )
                t = 1. / ( abs( theta ) + sqrt( 1. + theta ** 2 ) )
                if ( theta .lt. 0. ) t = -t
              endif ! abs(h)
              c     = 1. / sqrt( 1. + t**2 )
              s     = t * c
              tau   = s / ( 1. + c )
              h     = t * a( ip, iq )
              z( ip ) = z( ip ) - h
              z( iq ) = z( iq ) + h
              d( ip ) = d( ip ) - h
              d( iq ) = d( iq ) + h
              a( ip, iq ) = 0.
              do j = 1, ip - 1
                g = a( j, ip )
                h = a( j, iq )
                a( j, ip ) = g - s * ( h + g * tau )
                a( j, iq ) = h + s * ( g - h * tau )
              enddo ! j
              do j = ip + 1, iq - 1
                g = a( ip,  j )
                h = a(  j, iq )
                a( ip,  j ) = g - s * ( h + g * tau )
                a(  j, iq ) = h + s * ( g - h * tau )
              enddo ! j
              do j = iq + 1, nn
                g = a( ip, j )
                h = a( iq, j )
                a( ip, j ) = g - s * ( h + g * tau )
                a( iq, j ) = h + s * ( g - h * tau )
              enddo ! j
              do j = 1, nn
                g = v( j, ip )
                h = v( j, iq )
                v( j, ip ) = g - s * ( h + g * tau )
                v( j, iq ) = h + s * ( g - h * tau )
              enddo ! j
              m = m + 1
            endif ! i
          enddo ! iq
        enddo ! ip
        do ip = 1, nn
          b( ip ) = b( ip ) + z( ip )
          d( ip ) = b( ip )
          z( ip ) = 0.
        enddo ! ip
      enddo ! i

      stop 'Jacobi: too many iterations'

      end ! Jacobi
