B    = b*E
     a    = apply(B,1,sum)
     a    = max(a) + 1
     a    = rep(a,nrow(E))
     Sinv = B + diag(a)




draw = function(E){
     par(pty="s")
     n = nrow(E)
     angle = seq(0,2*pi,length=n+1)
     angle = angle[1:n]
     x = cos(angle)
     y = sin(angle)
     plot(x,y,pch=20,lwd=3,xlab="",ylab="",
          xaxt="n",yaxt="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
     for(i in 1:n){
          text(1.2*cos(angle[i]),1.2*sin(angle[i]),paste("X",i),lwd=3,font=2)
          }
     for(i in 1:(n-1)){
          for(j in (i+1):n){
               if(E[i,j] == 1)lines(c(x[i],x[j]),c(y[i],y[j]),lwd=3)
               }
          }
     return(NULL)
     }



     e = eigen(A)
     V = e$vectors
     s = V %*% sqrt(diag(e$values)) %*% t(V)