Dear KS,

On a linux workstation the ratio of runtimes is around 2.3 in favor of
MUMPS, so the problem seems to pertain to your installation or to the
windows installations in general. Can anyone using Kwant on windows
try the attached script? (I reduced the number of points so that it
doesn't take 1/2 hour and removed some parts not crucial to the
benchmark)

Best,
Anton

On Tue, May 16, 2017 at 7:18 PM, Prof. CHAN Kwok Sum
<apksc...@cityu.edu.hk> wrote:
> Dear Anton,
> To confirm my observation, I remove the non-essential part of the script to 
> test the mumps solver and default solver. The system is a zigzag graphene 
> ribbon, with 200X200 unit cells, there are 4 atoms in one unit cell so there 
> are 1.6x10^5 atoms. The OS is window 10, kwant is 1.2.2 and python is 3.4.4. 
> the time for calculating 20 smatrix is around 30 minutes for both default 
> solver and mumps. Please find attached the script and the runtime. The first 
> runtime is for mumps. The second runtime is for default. The kwant was 
> installed with the package prepared by Gohlke.
> For the sparse solver the time is longer, 40 minutes but there is no umfpack 
> installed in the window system. So I do not show the runtime as it may not be 
> accurate.
> KS
>
> On 16/5/2017, 5:51 PM, "anton.akhme...@gmail.com on behalf of Anton Akhmerov" 
> <anton.akhme...@gmail.com on behalf of anton.akhmerov...@gmail.com> wrote:
>
>     Dear KS,
>
>     > I have compared the performance of default, mumps solvers in windows and
>     > Ubuntu machines. Their performance is very similar without significant 
> gain
>     > when using mumps. Do I miss anything in the installation? For example, 
> not
>     > installing the mumps. Similar performance is found for the sparse 
> solver.
>     > The system I considered has 10^5 to 6X10^5 sites and machine has a quad 
> core
>     > i7 cpu with 16G RAM. My installation procedures are those found in the
>     > website for pre complled packages.
>
>     This is contrary to our experience. Can you please provide more details:
>     * How did you install Kwant in both cases?
>     * What specific script did you run?
>     * What were the execution times that you observed?
>
>     Thanks,
>     Anton
>
>     > KS
>     >
>     >
>     >
>     > Disclaimer: This email (including any attachments) is for the use of the
>     > intended recipient only and may contain confidential information and/or
>     > copyright material. If you are not the intended recipient, please 
> notify the
>     > sender immediately and delete this email and all copies from your 
> system.
>     > Any unauthorized use, disclosure, reproduction, copying, distribution, 
> or
>     > other form of unauthorized dissemination of the contents is expressly
>     > prohibited.
>
>
>
>
> Disclaimer: This email (including any attachments) is for the use of the 
> intended recipient only and may contain confidential information and/or 
> copyright material. If you are not the intended recipient, please notify the 
> sender immediately and delete this email and all copies from your system. Any 
> unauthorized use, disclosure, reproduction, copying, distribution, or other 
> form of unauthorized dissemination of the contents is expressly prohibited.
import os
from time import time
from math import pi, sqrt, tanh
from cmath import exp


import kwant
import kwant.solvers.sparse as ksparse
import kwant.solvers.mumps as kmumps
import numpy as np


###setup directory
##curdir=os.getcwd()
##os.chdir(curdir+'/test')
print (ksparse.uses_umfpack)
#define function for magnetic field phase shift
def hopphase_x(sitei, sitej, flux):

    xi, yi = sitei.pos
    xj, yj = sitej.pos
    return exp(-0.5*1j*flux*(xi-xj)*(yi+yj))
def hopphase_y(sitei, sitej, flux):
     xi, yi = sitei.pos
     xj, yj = sitej.pos
     return exp(0.5*1j* flux * (xi + xj) * (yi - yj))

#HIDDEN_END_dwhx
print ('start',time())
sys = kwant.Builder()
sys2=kwant.Builder()
#HIDDEN_END_goiq

# Here, we are only working with square lattices
#HIDDEN_BEGIN_suwo
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
#a = 1
# lat = kwant.lattice.honeycomb(a)
# this lattice is not a honeycomb lattice. it is like a rectangle
# this lattice gives correct zigzag and armchair edges
graphene= kwant.lattice.general([(sqrt(3), 0), (0, 1)],
                           [(0,1/2), (1/(2*sqrt(3)), 0),(sqrt(3)/2, 0), (2/sqrt(3),1/2)])



def rect(pos):
    x, y= pos
    lx=L*sqrt(3);ly=W*1;
    return 0<=x<lx and 0<=y<ly

#Lstruct and Lshift are labels for file name indicating the structure and shift
Lstruct="_W40L40"
Lshift="no shift"
pp = PdfPages('ribbon'+Lstruct+Lshift+'.pdf')
t = 1.0
W = 200   # y dimension scatteribg region and horizontal lead
Vposit=0    #position of the horizontal lead
WY=200  # width of horizontal lead which can be different from W
L = 200   # x dimension scattering region
WL=200     #WL is width of vertical lead; it could be different from L
lpos=0   # position of vertical lead; the L and WL are different
lpos2=0
Lstub=0 #length of vertical stub which can have potential
WL1=2
WL2=7
flux=0.00
shiftcenter=0.00 #shift the center scattering region
shiftlead=-0.0
enpt=0.3 # shift the energy range.


a, b, c, d =graphene.sublattices

# Define the scattering region
sys2[graphene.shape(rect,(0.,0.))]=0.0*t
sys2[graphene.neighbors()]=-t
for i in range(L):
    wtmp=0 #wtmp define the width of the central ribbon
    if lpos+WL-1>=i>=lpos: wtmp=-Lstub
    for j in range(wtmp,W):
        # On-site Hamiltonian
#     if i<23 : 
#        sys[a(i, j)] = 0. * t
#        sys[b(i,j)]=0.*t
#        sys[c(i,j)]=0.*t
#        sys[d(i,j)]=0.*t
#      else :
        sys[a(i, j)] = (0.+i*0.0)* t+shiftcenter
        sys[b(i,j)]=(0.+i*0.0)*t+shiftcenter
        sys[c(i,j)]=(0.+i*0.0)*t+shiftcenter
        sys[d(i,j)]=(0.+i*0.0)*t+shiftcenter
#new hopping
        
        sys[a(i,j),b(i,j)]=-t*hopphase_x(a(i,j),b(i,j),flux)
        sys[b(i,j),c(i,j)]=-t*hopphase_x(b(i,j),c(i,j),flux)
        sys[c(i,j),d(i,j)]=-t*hopphase_x(c(i,j),d(i,j),flux)
 
        if j<0: icheck=lpos
        else : icheck=0
        if i>icheck:
#          (xa,ya)=a(i,j).pos
#          (xd,yd)=d(i-1,j).pos
#          sys[a(i,j), d(i-1,j)]=-t*exp(flux*1j*0.5*(xa-xd)*(ya+yd))
          sys[a(i,j), d(i-1,j)]=-t*hopphase_x(a(i,j) ,d(i-1,j),flux)        
        if lpos+WL>=i>=lpos: jcheck=-Lstub
        else : jcheck=0
        if j>jcheck:
##           (xa,ya)=a(i,j-1).pos
##           (xb,yb)=b(i,j).pos
##           sys[b(i,j), a(i,j-1)]=-t*exp(flux*1j*0.5*(xb-xa)*(ya+yb))
           sys[b(i,j), a(i,j-1)]=-t*hopphase_x(b(i,j),a(i,j-1),flux)
##           (xc,yc)=c(i,j-1).pos
##           (xd,yd)=d(i,j).pos            
##           sys[c(i,j),d(i,j-1)]=-t*exp(flux*1j*0.5*(xc-xd)*(yc+yd))
           sys[c(i,j), d(i,j-1)]=-t*hopphase_x(c(i,j),d(i,j-1),flux)      



# vertical lead
sym2 = kwant.TranslationalSymmetry(graphene.vec((0,-1)))
lead2=kwant.Builder(sym2)
for i in range(lpos, WL+lpos):
    lead2[a(i,0)]=0.*t
    lead2[b(i,0)]=0.*t
    lead2[c(i,0)]=0.*t
    lead2[d(i,0)]=0.*t
#    lead2[graphene.neighbors()]=-t
    lead2[a(i,0),b(i,0)]=-t*hopphase_y(a(i,0),b(i,0),flux)
    lead2[b(i,0),c(i,0)]=-t*hopphase_y(b(i,0),c(i,0),flux)
    lead2[c(i,0),d(i,0)]=-t*hopphase_y(c(i,0),d(i,0),flux)
    lead2[b(i,0),a(i,-1)]=-t*hopphase_y(b(i,0),a(i,-1),flux)
    lead2[c(i,0),d(i,-1)]=-t*hopphase_y(c(i,0),d(i,-1),flux)
 
    if i>lpos:
      lead2[a(i,0),d(i-1,0)]=-t*hopphase_y(a(i,0),d(i-1,0),flux)
 
sym3=kwant.TranslationalSymmetry(graphene.vec((0,1)))
lead3=kwant.Builder(sym3)
for i in range(lpos2, WL+lpos2):
    lead3[a(i,0)]=0.0*t
    lead3[b(i,0)]=0*t
    lead3[c(i,0)]=0.0*t
    lead3[d(i,0)]=0*t
#    lead1[graphene.neighbors()]=-t
    lead3[a(i,0),b(i,0)]=-t*hopphase_y(a(i,0),b(i,0),flux)
    lead3[b(i,0),c(i,0)]=-t*hopphase_y(b(i,0),c(i,0),flux)
    lead3[c(i,0),d(i,0)]=-t*hopphase_y(c(i,0),d(i,0),flux)
    lead3[b(i,0),a(i,-1)]=-t*hopphase_y(b(i,0),a(i,-1),flux)
    lead3[c(i,0),d(i,-1)]=-t*hopphase_y(c(i,0),d(i,-1),flux)
    if i>lpos2:
      lead3[a(i,0),d(i-1,0)]=-t*hopphase_y(a(i,0),d(i-1,0),flux)
       
#attach the lead here

sys.attach_lead(lead2)
sys.attach_lead(lead3)


# plot band structure of lead2
lead2plot = lead2.finalized()

kwant.plotter.bands(lead2plot, show=False)

fsys = sys.finalized()

energies = []
data = []  #data array for lead 2
data_chck=[]


t = -time()

for ie in range(-650, 650, 250):
    energy = ie * 0.0005+0.00001+enpt
    print (ie, time() + t)
#    print energy

    # compute the scattering matrix at a given energy
    smatrix = kwant.smatrix(fsys, energy)


print ('MUMPS: ', time() + t )

t = -time()

for ie in range(-650, 650, 250):
    energy = ie * 0.0005+0.00001+enpt
    print (ie, time() + t)
#    print energy

    # compute the scattering matrix at a given energy
    smatrix = kwant.solvers.sparse.smatrix(fsys, energy)


print ('SuperLU: ', time() + t )

Reply via email to