hmat_div.tcell

(%ifndef* NF-TYPE
  (%defconstant NF-TYPE GCC)) ; one of (GCC LW-SC CL-SC XCC XCCCL) 
(%include "rule/tcell-setrule.sh")

(c-exp "#include <stdio.h>")
(c-exp "#include <stdlib.h>")
(c-exp "#include <math.h>")
(c-exp "#include <time.h>")
(c-exp "#include <sys/time.h>")
(c-exp "#include <unistd.h>")
(c-exp "#include \"../data/bem_file.h\"")
(%defconstant PN 10000)             ;;TS
(%defconstant PL 10)                ;;TN
(%defconstant PARA_LEVEL 20)
(%defconstant CHUNK_SIZE 128)      ;;C
(%defconstant Check-interval 100)
(%defconstant MB 25000000)

(%defconstant TCELL-CHECK-REQ 1)

;; ; Spherecube
;;(%defconstant INPUT_DEFAULT "../data_pro1804/input_100ts.txt_cb_1.5_10_10_10.bin")
;; ; Spherepyramid
;;(%defconstant INPUT_DEFAULT "../data_pro1804/input_100ts.txt_pb_1.5_14.bin")
;; ; Sphere
(%defconstant INPUT_DEFAULT "../data_pro1804/input_50ms.bin")
;; ; humanoids
;(%defconstant INPUT_DEFAULT "../data_pro1804/input_human_1x1.txt_cb_0.3_50_100_1.bin")
;; ; small dataset for debug, N=1000
;(%defconstant INPUT_DEFAULT "../data_pro1804/input_100ts.txt.bin")

(%defconstant NWORKER-MAX 540)
(extern-decl num-thrs unsigned-int)

;; *********define cluster************
(deftype cluster (struct cluster))
(def (struct cluster)
  (decl ndim int)
  (decl nstrt int)
  (decl nsize int)
  (decl ndpth int)
  (decl nnson int)
  (decl nmbr int)
  (decl ndscd int)
  (decl bmin (ptr double))
  (decl bmax (ptr double))
  (decl zwdth double)
  (decl pc-sons (ptr (ptr cluster))))

(deftype leafmtx (struct leafmtx))
(def (struct leafmtx)
  (decl ltmtx int)
  (decl kt int)
  (decl nstrtl int)
  (decl ndl int)
  (decl nstrtt int)
  (decl ndt int)
  (decl a1 (ptr double))
  (decl a2 (ptr double)))

(deftype leafmtxp (struct leafmtxp))
(def (struct leafmtxp)
  (decl nlf int)
  (decl nlfkt int))

(def worker-data
  (def n int)
  (def leafnode-list (ptr leafmtx))
  (def lnmtx (array int 3))
  (def nlf int))

;; Pointer to WDATA of i-th worker
(def Wdata-list (array (ptr (struct worker-data)) NWORKER-MAX))

(def worker-init
  (= (aref Wdata-list WORKER-ID) (ptr WDATA))
  (= WDATA.n 200000000)
  ;; Do malloc() after CT construction
  ;;(= WDATA.leafnode-list (cast (ptr leafmtx) (csym::malloc (* WDATA.n (sizeof leafmtx)))))
  (= (aref WDATA.lnmtx 0) 0)
  (= (aref WDATA.lnmtx 1) 0)
  (= (aref WDATA.lnmtx 2) 0)
  (= WDATA.nlf 0))

(def (task create-leafmtx)
  (def st-cltl (ptr cluster) :in)
  (def st-cltt (ptr cluster) :in)
  (def i1 int :in)
  (def i2 int :in))

(def (task start-up)
  (def r int))

(def (task build-cluster)
  ;;(def st-clt (ptr cluster) :in)
  (def zgmid (ptr (array double 3)) :in)
  (def tempzgmid (ptr (array double 3)) :in)
  (def ndpth int :in)
  (def ndscd int :in)
  (def nsrt int :in)
  (def nd int :in)
  (def md int :in)
  (def nclst int :in)
  (def r (ptr cluster))
)

(def (task minmax)
  (def i1 int :in)
  (def i2 int :in)
  (def nsrt int :in)
  (def zgmid (ptr (array double 3)) :in)
  (def ndim int :in)
  (def ndpth int :in)
  (def minx double)
  (def maxx double)
  (def miny double)
  (def maxy double)
  (def minz double)
  (def maxz double))

(def (task countNum)
  (def i1 int :in)
  (def i2 int :in)
  (def zgmid (ptr (array double 3)) :in)
  (def nd int :in)
  (def ncut int :in)
  (def zlmid double :in)
  (def gn int :in)
  (def lessNum (ptr int) :in)
  (def moreNum (ptr int) :in)
  (def ndpth int :in)
  (def nsrt int :in)
)

(def (task setStart)
  (def i1 int :in)
  (def i2 int :in)
  (def nd int :in)
  (def ncut int :in)
  (def gn int :in)
  (def nsrt int :in)
  (def zlmid double :in)
  (def zgmid (ptr (array double 3)) :in)
  (def tempzgmid (ptr (array double 3)) :in)
  (def lessNum (ptr int) :in)
  (def moreNum (ptr int) :in)
  (def lessStart (ptr int) :in)
  (def moreStart (ptr int) :in)
  (def ndpth int :in))

(def (task ssum)
  (def i1 int :in)
  (def i2 int :in)
  (def lessNum (ptr int) :in)
  (def ndpth int :in)
  (def sum int)
)

(decl (start-up) (wfn int))
(decl (supermatrix-construction-cog-leafmtrx st-leafmtxp gmid param lod lnmtx nofc nffc ndim)
    (wfn void (ptr leafmtxp) (ptr (array double 3)) (array double) (ptr int) (ptr int) int int int))
(decl (csym::qsort-col-leafmtx st-leafmtx first last)
    (fn void (ptr leafmtx) int int))
(decl (csym::qsort-row-leafmtx st-leafmtx first last)
    (fn void (ptr leafmtx) int int))
(decl (csym::med3 nl nr nlr2)
    (fn int int int int))
(decl (create-leafmtx stc-cltl st-cltt i1 i2)
    (wfn void (ptr cluster) (ptr cluster) int int))
(decl (csym::dist-2cluster st-cltl st-cltt)
    (fn double (ptr cluster) (ptr cluster)))
(decl (csym::count-lntmx st-cltl st-cltt param lnmtx nffc)
    (fn void (ptr cluster) (ptr cluster) (array double) (ptr int) int))
(decl (csym::cal-bndbox-cog st-clt zgmid nofc)
    (fn void (ptr cluster) (ptr (array double 3)) int))
(decl (csym::set-bndbox-cog st-clt zgmid nofc)
    (fn void (ptr cluster) (ptr (array double 3)) int))
(decl (csym::create-cluster nmbr ndpth nstrt nsize ndim nson)
    (fn (ptr cluster) int int int int int int))
(decl (csym::free-st-clt st-clt)
    (fn void (ptr cluster)))
(decl (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt nd md nclst)
    (wfn (ptr cluster) (ptr (array double 3)) (ptr (array double 3)) int int int int int int))
(decl (csym::get-wall-time) (fn double))
(decl (csym::get-cpu-time) (fn double))
(decl (csym::checkClusterTree f st-clt) (fn void (ptr FILE) (ptr cluster)))
(decl (minmax i1 i2 zgmid tsk ndpth nsrt) (wfn void int int (ptr (array double 3)) (ptr (struct minmax)) int int))
(decl (countNum i1 i2 lessNum moreNum zgmid nd ncut glmid gn ndpth nsrt) (wfn void int int (ptr int) (ptr int) (ptr (array double 3)) int int double int int int))
(decl (setStart i1 i2 ncut gn nd zlmid lessNum moreNum lessStart moreStart zgmid tempzgmid ndpth nsrt)
    (wfn void int int int int int double (ptr int) (ptr int) (ptr int) (ptr int) (ptr (array double 3)) (ptr (array double 3)) int int))
(decl (ssum i1 i2 ndpth lessNum) (wfn int int int int (ptr int)))

(decl depth-max int)
(decl count-node int)
(decl param (array double 100))
(decl nffc int)
(decl nff (array int 100))
(decl ndim int)

(def (task-body create-leafmtx)
  (create_leafmtx this.st-cltl this.st-cltt this.i1 this.i2))

(def (task-body start-up)
  (= this.r (start-up))
)

(def (task-body build-cluster)
  (= this.r (create-ctree-ssgeom this.zgmid this.tempzgmid this.ndpth 
				 this.ndscd this.nsrt this.nd this.md this.nclst)))

(def (task-body minmax)
  (minmax this.i1 this.i2 this.zgmid (ptr this) this.ndpth this.nsrt))

(def (task-body countNum)
  ;;(csym::fprintf stderr "countNum: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID this.i1 this.i2 this.ndpth)
  (countNum this.i1 this.i2 this.lessNum this.moreNum this.zgmid
	    this.nd this.ncut this.zlmid this.gn this.ndpth this.nsrt))

(def (task-body setStart)
  ;;(csym::fprintf stderr "setStart: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID this.i1 this.i2 this.ndpth)
  (setStart this.i1 this.i2 this.ncut this.gn this.nd this.zlmid this.lessNum this.moreNum
	    this.lessStart this.moreStart this.zgmid this.tempzgmid this.ndpth this.nsrt))

(def (task-body ssum)
  (= this.sum (ssum this.i1 this.i2 this.ndpth this.lessNum))
)

(def (start-up) (wfn int)
  (decl fname (ptr char))
  (decl file (ptr FILE))
  (def countOfNode int 0)
  (def count int 0)
  (decl i int)
  (def bi (struct bem-input))
  (decl coordOfNode (ptr (array double 3)))
  (decl coordOfFace (ptr (array double 3)))
  (= fname INPUT-DEFAULT)
  (= file (fopen fname "r"))
  (if (== file NULL)
      (begin
        (csym::fprintf stderr "Error: Unable to input file '%s'!~%" fname)
        (csym::exit 99))
    (begin
      (csym::fprintf stderr "Input file: %s~%" INPUT_DEFAULT)
      (if (== (csym::read-bem-input file (ptr bi) BI-AUTO) -1)
	  (begin
	    (csym::fprintf stderr "Bem input file read error!~%")
	    (csym::exit 99)))
      (csym::print-bem-input stderr (ptr bi) BI-PRETTY)
      (= countOfNode bi.nNode)
      (= coordOfNode bi.coordOfNode)
      (= count bi.nFace)
      (= coordOfFace bi.coordOfFace)
      ))
  ;(csym::fprintf stderr "1~%")
  (csym::fclose file)
  (csym::free coordOfNode)
  (for ((= i 0) (< i 100) (inc i))
    (= (aref param i) 0.0)
    (= (aref nff i) 0))
  (= (aref param 21) 10.0)
  (= (aref param 31) 1.1)
  (= (aref param 41) 15.0)
  (= (aref param 51) 2.0)
  (decl st-leafmtxp (ptr leafmtxp))
  (decl lod (ptr int))
  (decl lnmtx (ptr int))
  (def nofc int count)
  (= nffc 1)
  ;(def ndim int 3)
  (= ndim 3)
  (= lnmtx (cast (ptr int) (csym::malloc (* 3 (sizeof int)))))
  (for ((= i 0) (< i 3) (inc i))
    (= (aref lnmtx i) 0))
  (= st-leafmtxp (cast (ptr leafmtxp) (csym::malloc (sizeof leafmtxp))))
  (= lod (cast (ptr int) (csym::malloc (* nofc (sizeof int)))))
  (for ((= i 0) (< i nofc) (inc i))
    (= (aref lod i) 0))
  ;(csym::fprintf stderr "2~%")
  (supermatrix-construction-cog-leafmtrx st-leafmtxp coordOfFace param lod lnmtx
                                         nofc nffc ndim)
  (return 0))

(def (supermatrix-construction-cog-leafmtrx st-leafmtxp gmid param lod lnmtx nofc nffc ndim)
    (wfn void (ptr leafmtxp) (ptr (array double 3)) (array double) (ptr int) (ptr int) int int int)
  ;;(def st-clt (ptr cluster) (cast (ptr cluster) (csym::malloc (sizeof cluster))))
  (def st-clt (ptr cluster))
  (decl i int)
  (decl nfl int)
  (decl nflkt int)
  (decl ip int) 
  (decl il int)
  (decl ig int)
  (def nd int (* nofc nffc))
  (decl lodfc (ptr int))
  (decl st-leafmtx (ptr leafmtx))
  (= lodfc (cast (ptr int) (csym::malloc (* nofc (sizeof int)))))
  (for ((= il 0) (< il nofc) (inc il))
    (= (aref lodfc il) il))
  (def nsrt int 0)
  (def ndf int nofc)
  (def nclst int 0)
  (def ndpth int 0)
  (def ndscd int 0)
  (= depth-max 0)
  (= count-node 0)
  (decl start double)
  (decl end double)
  (decl spent double)

  (decl tempgmid (ptr (array double 3)))
  ;(csym::fprintf stderr "3~%")
  (= tempgmid (cast (ptr (array double 3))
		(csym::malloc (* (* nofc 3) (sizeof double)))))
  ;(csym::fprintf stderr "4~%")

  ;(csym::fprintf stderr "test!~%")
  (= start (csym::get-wall-time))
  (= st-clt
     (create-ctree-ssgeom gmid tempgmid ndpth ndscd nsrt
                          ndf nofc nclst))
  (= end (csym::get-wall-time))
  (= spent (- end start))
  (csym::fprintf stderr "cluster tree time spent:%.10lf~%" spent)
  ;;(csym::set-bndbox-cog st-clt gmid nofc)

  (decl f (ptr FILE))
  (= f (csym::fopen "hmat_all.txt" "w"))
  ;;(csym::sleep 10)
  (checkClusterTree f st-clt)
  ;;(csym::sleep 10)
  (csym::fclose f)

  (= ndpth 0)
  (= start (csym::get-wall-time))
  (csym::count-lntmx st-clt st-clt param lnmtx nffc)
  (= end (csym::get-wall-time))
  (= spent (- end start))
  (csym::fprintf stderr "count time:%.10lf~%" spent)
  (= (fref (mref st-leafmtxp) nlfkt) (aref lnmtx 0))
  (def nlf int (+ (aref lnmtx 0) (aref lnmtx 1)))
  (= st-leafmtx (cast (ptr leafmtx) (csym::malloc (* nlf (sizeof leafmtx)))))
  (= (fref (mref st-leafmtxp) nlf) nlf)
  (csym::fprintf stderr "nlf: %ld~%" nlf)
  
  (def nlf-for-each-worker int (/ (* (cast long nlf) 100) num-thrs))
  (csym::fprintf stderr "nlf-for-each-worker: %ld num-thrs: %ld ~%" nlf-for-each-worker num-thrs)
  (begin
    (def sz size-t (* nlf-for-each-worker (cast long num-thrs)))
    (for ((= i 0) (< i num-thrs) (inc i))
      (= (fref (aref Wdata-list i) -> leafnode-list)
         (cast (ptr leafmtx) (csym::malloc sz)))
      (if (not (fref (aref Wdata-list i) -> leafnode-list))
          (begin
            (csym::fprintf stderr "malloc temp-leafmtx for worker %d error!~%" i)
            (csym::exit 99)))))
  (= nlf 0)
  (decl all int)
  (= all (* (fref (mref st-clt) nnson) (fref (mref st-clt) nnson)))
  (csym::fprintf stderr "all:%ld~%" all)
  (= start (csym::get-wall-time))
  (create-leafmtx st-clt st-clt 0 all)
  (= end (csym::get-wall-time))
  (= spent (- end start))

  (def nncon int 0)
  (for ((= i 0) (< i 100) (inc i))
    (= nncon (+ (aref nff i) nncon)))
  ;(csym::fprintf stderr "nlf:%ld~%" WDATA.nlf)
  (begin
    (def nlf-sum int 0)
    (for ((= i 0) (< i num-thrs) (inc i))
      (def nlf-i int (fref (aref Wdata-list i) -> nlf))
      (csym::fprintf stderr "nlf(%d) = %d~%" i nlf-i)
      (if (> nlf-i nlf-for-each-worker)
          (begin
            (csym::fprintf stderr "nlf of counted by worker %d is larger than nlf-for-each-worker!~%" i)
            (csym::exit 99)))
      (+= nlf-sum nlf-i))
    (csym::fprintf stderr "nlf-sum = %d~%" nlf-sum))
  
  ;;(csym::fprintf stderr "nlf:%ld~%" nncon)
  (csym::fprintf stderr "block cluster tree time spent:%.10lf~%" spent)
  (csym::fprintf stderr "depth_max:%ld  count_node:%ld~%" depth-max count-node))

(def (csym::qsort-row-leafmtx st-leafmtx first last)
    (fn void (ptr leafmtx) int int)
  (decl i int) 
  (decl j int)
  (decl pivot int)
  (decl st-www leafmtx)
  (if (< first last)
      (begin 
        (= pivot first)
        (= i first)
        (= j last)
        (while (< i j)
          (while (and (<= (fref (aref st-leafmtx i) nstrtl)
                          (fref (aref st-leafmtx pivot) nstrtl))
                      (< i last))
            (inc i))
          (while (> (fref (aref st-leafmtx j) nstrtl)
                    (fref (aref st-leafmtx pivot) nstrtl))
            (dec j))
          (if (< i j)
              (begin
                (= st-www (aref st-leafmtx i))
                (= (aref st-leafmtx i) (aref st-leafmtx j))
                (= (aref st-leafmtx j) st-www))))
        (= st-www (aref st-leafmtx pivot))
        (= (aref st-leafmtx pivot) (aref st-leafmtx j))
        (= (aref st-leafmtx j) st-www)
        (csym::qsort-row-leafmtx st-leafmtx first (- j 1))
        (csym::qsort-row-leafmtx st-leafmtx (+ j 1) last))))

(def (csym::qsort-col-leafmtx st-leafmtx first last) (fn void (ptr leafmtx) int int)
  (decl i int)
  (decl j int)
  (decl pivot int)
  (decl st-www leafmtx)
  (if (< first last)
      (begin 
        (= pivot first)
        (= i first)
        (= j last)
        (while (< i j)
          (while
              (and (<= (fref (aref st-leafmtx i) nstrtt)
                       (fref (aref st-leafmtx pivot) nstrtt))
                   (< i last))
            (inc i))
          (while
              (> (fref (aref st-leafmtx j) nstrtt) (fref (aref st-leafmtx pivot) nstrtt))
        (dec j))
       (if (< i j)
           (begin (= st-www (aref st-leafmtx i))
            (= (aref st-leafmtx i) (aref st-leafmtx j))
            (= (aref st-leafmtx j) st-www))))
      (= st-www (aref st-leafmtx pivot))
      (= (aref st-leafmtx pivot) (aref st-leafmtx j))
      (= (aref st-leafmtx j) st-www)
      (csym::qsort-col-leafmtx st-leafmtx first (- j 1))
      (csym::qsort-col-leafmtx st-leafmtx (+ j 1) last))))

(def (csym::med3 nl nr nlr2) (fn int int int int)
  (decl med3 int)
  (if (< nl nr)
      (begin
        (if (< nr nlr2)
            (begin (= med3 nr))
            (if (< nlr2 nl)
                (begin (= med3 nl))
                (begin (= med3 nlr2)))))
      (begin
        (if (< nlr2 nr)
            (begin (= med3 nr))
            (if (< nl nlr2) 
                (begin (= med3 nl))
                (begin (= med3 nlr2))))))
  (return med3))

(def (create-leafmtx st-cltl st-cltt i1 i2)
    (wfn void (ptr cluster) (ptr cluster) int int)
  (def ndl int (* (fref (mref st-cltl) nsize) nffc))
  (def ndt int (* (fref (mref st-cltt) nsize) nffc))
  (def nstrtl int (fref (mref st-cltl) nstrt))
  (def nstrtt int (fref (mref st-cltt) nstrt))
  (def nnsonl int (fref (mref st-cltl) nnson))
  (def nnsont int (fref (mref st-cltt) nnson))
  (def ndpth int (fref (mref st-cltt) ndpth))
  (decl ia int)
  (def nleaf double (aref param 41))
  (def zeta double (aref param 51))
  (def zdistlt double (csym::dist-2cluster st-cltl st-cltt))
  (if (and (or (<= (* (fref (mref st-cltl) zwdth) zeta) zdistlt)
               (<= (* (fref (mref st-cltt) zwdth) zeta) zdistlt))
           (and (>= ndl nleaf) (>= ndt nleaf)))
      (begin
	;;(csym::fprintf stderr "type1 nlf:%ld~%" WDATA.nlf)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) nstrtl) nstrtl)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) ndl) ndl)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) nstrtt) nstrtt)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) ndt) ndt)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) kt) 0)
        (= (fref (aref WDATA.leafnode-list WDATA.nlf) ltmtx) 1)
        (= WDATA.nlf (+ WDATA.nlf 1))
	;(csym::fprintf stderr "%ld end~%" WDATA.nlf)
	;(csym::fprintf stderr "type1 end~%")
	;;(= (aref nff WORKER-ID) (+ (aref nff WORKER-ID) 1))
	#+comment(csym::fprintf stderr "type1 end~%"))
      (begin
        (if (or (== nnsonl 0) (== nnsont 0) (<= ndl nleaf) (<= ndt nleaf))
            (begin
	      ;(csym::fprintf stderr "type2 nlf:%ld~%" WDATA.nlf)
	      (= (fref (aref WDATA.leafnode-list WDATA.nlf) nstrtl) nstrtl)
              (= (fref (aref WDATA.leafnode-list WDATA.nlf) ndl) ndl)
              (= (fref (aref WDATA.leafnode-list WDATA.nlf) nstrtt) nstrtt)
              (= (fref (aref WDATA.leafnode-list WDATA.nlf) ndt) ndt)
              (= (fref (aref WDATA.leafnode-list WDATA.nlf) ltmtx) 2)
              (= WDATA.nlf (+ WDATA.nlf 1))
	      ;(csym::fprintf stderr "%ld end~%" WDATA.nlf)
	      ;;(= (aref nff WORKER-ID) (+ (aref nff WORKER-ID) 1))
	      #|(csym::fprintf stderr "type2 end~%")|#)
            (begin
	      (if (< ndpth PARA_LEVEL)
		  (begin
		    ;(csym::fprintf stderr "type3 end~%")
		    (do-many for ia from i1 to i2
		      (decl il int) (decl it int)
		      (= il (/ ia nnsont))
		      (= it (% ia nnsont))
		      (decl al int)
		      (= al (* (fref (aref (fref st-cltl -> pc-sons) il) -> nnson) (fref (aref (fref st-cltt -> pc-sons) it) -> nnson)))
		      (create-leafmtx (aref (fref (mref st-cltl) pc-sons) il) (aref (fref (mref st-cltt) pc-sons) it) 0 al)
		      (handles create-leafmtx
			(:put from j1 to j2
			      (= this.st-cltl st-cltl)
			      (= this.st-cltt st-cltt)
			      (= this.i1 j1)
			      (= this.i2 j2))
			(:get)))
		    )
		(begin
		  ;(csym::fprintf stderr "type4 end~%")
		  (for ((= ia i1) (< ia i2) (inc ia))
		    (decl il int) (decl it int)
		    (= il (/ ia nnsont))
		    (= it (% ia nnsont))
		    (decl al int)
		    (= al (* (fref (aref (fref st-cltl -> pc-sons) il) -> nnson) (fref (aref (fref st-cltt -> pc-sons) it) -> nnson)))
		    (create-leafmtx (aref (fref (mref st-cltl) pc-sons) il) (aref (fref (mref st-cltt) pc-sons) it) 0 al)
		    )
		  )
		)
	      )))))

(def (csym::dist-2cluster st-cltl st-cltt) (fn double (ptr cluster) (ptr cluster))
 (def zs double 0.0) (decl id int)
 (for ((= id 0) (< id (fref (mref st-cltl) ndim)) (inc id))
   (if (< (aref (fref (mref st-cltl) bmax) id) (aref (fref (mref st-cltt) bmin) id))
       (begin
         (= zs
            (+ zs
               (* (- (aref (fref (mref st-cltt) bmin) id)
                     (aref (fref (mref st-cltl) bmax) id))
                  (- (aref (fref (mref st-cltt) bmin) id)
                     (aref (fref (mref st-cltl) bmax) id))))))
       (if (< (aref (fref (mref st-cltt) bmax) id)
              (aref (fref (mref st-cltl) bmin) id))
           (begin
             (= zs
                (+ zs
                   (* (- (aref (fref (mref st-cltl) bmin) id)
                         (aref (fref (mref st-cltt) bmax) id))
                      (- (aref (fref (mref st-cltl) bmin) id)
                         (aref (fref (mref st-cltt) bmax) id)))))))))
 (return (sqrt zs)))

(def (csym::count-lntmx st-cltl st-cltt param lnmtx nffc)
    (fn void (ptr cluster) (ptr cluster) (array double) (ptr int) int)
  (decl il int)
  (decl it int)
  (def ndl int (* (fref (mref st-cltl) nsize) nffc))
  (def ndt int (* (fref (mref st-cltt) nsize) nffc))
  (def nstrtl int (fref (mref st-cltl) nstrt))
  (def nstrtt int (fref (mref st-cltt) nstrt))
  (def nnsonl int (fref (mref st-cltl) nnson))
  (def nnsont int (fref (mref st-cltt) nnson))
  (def nleaf double (aref param 41))
  (def zeta double (aref param 51))
  (def zdistlt double (dist-2cluster st-cltl st-cltt))
  (if (and (or (<= (* (fref (mref st-cltl) zwdth) zeta) zdistlt)
               (<= (* (fref (mref st-cltt) zwdth) zeta) zdistlt))
           (and (>= ndl nleaf) (>= ndt nleaf)))
      (begin
        (= (aref lnmtx 0) (+ (aref lnmtx 0) 1)))
      (begin
        (if (or (== nnsonl 0) (== nnsont 0) (<= ndl nleaf) (<= ndt nleaf))
            (begin
              (= (aref lnmtx 1) (+ (aref lnmtx 1) 1)))
            (begin
              (= (aref lnmtx 2) (+ (aref lnmtx 2) 1))
              (for ((= il 0) (< il nnsonl) (inc il))
                (for ((= it 0) (< it nnsont) (inc it))
                  (csym::count-lntmx (aref (fref (mref st-cltl) pc-sons) il)
                               (aref (fref (mref st-cltt) pc-sons) it) param lnmtx nffc))))))))

(def (csym::cal-bndbox-cog st-clt zgmid nofc)
    (fn void (ptr cluster) (ptr (array double 3)) int)
  (def ndim int (fref (mref st-clt) ndim))
  (decl id int)
  (decl il int)
  (= (fref (mref st-clt) bmin) (cast (ptr double) (csym::malloc (* 3 (sizeof double)))))
  (= (fref (mref st-clt) bmax) (cast (ptr double) (csym::malloc (* 3 (sizeof double)))))
  (def zeps double 1.0e-5)
  (if (> (fref (mref st-clt) nnson) 0)
      (begin
        (for ((= id 0) (< id ndim) (inc id))
          (= (aref (fref (mref st-clt) bmin) id)
             (aref (fref (mref (aref (fref (mref st-clt) pc-sons) 0)) bmin) id))
          (= (aref (fref (mref st-clt) bmax) id)
             (aref (fref (mref (aref (fref (mref st-clt) pc-sons) 0)) bmax) id)))
        (for ((= il 1) (< il (fref (mref st-clt) nnson)) (inc il))
          (for ((= id 0) (< id ndim) (inc id))
            (if (< (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmin) id)
                   (aref (fref (mref st-clt) bmin) id))
                (begin
                  (= (aref (fref (mref st-clt) bmin) id)
                     (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmax) id))))
            (if (< (aref (fref (mref st-clt) bmax) id)
                   (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmax) id))
                (begin
                  (= (aref (fref (mref st-clt) bmax) id)
                     (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmax) id)))))))
      (begin
        (for ((= id 0) (< id ndim) (inc id))
          (= (aref (fref (mref st-clt) bmin) id) (aref (aref zgmid 0) id))
          (= (aref (fref (mref st-clt) bmax) id) (aref (aref zgmid 0) id)))
        (for ((= id 0) (< id ndim) (inc id))
          (for ((= il 1) (< il (fref (mref st-clt) nsize)) (inc il))
            (if (< (aref (aref zgmid il) id)
                   (aref (fref (mref st-clt) bmin) id))
                (begin
                  (= (aref (fref (mref st-clt) bmin) id)
                     (aref (aref zgmid il) id))))
            (if (< (aref (fref (mref st-clt) bmax) id)
                   (aref (aref zgmid il) id))
                (begin
                  (= (aref (fref (mref st-clt) bmax) id)
                     (aref (aref zgmid il) id))))))))
  (def zwdth double
    (* (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))
       (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))))
  (for ((= id 1) (< id ndim) (inc id))
    (= zwdth
       (+ zwdth
          (* (- (aref (fref (mref st-clt) bmax) id)
                (aref (fref (mref st-clt) bmin) id))
             (- (aref (fref (mref st-clt) bmax) id)
                (aref (fref (mref st-clt) bmin) id))))))
  (= zwdth (csym::sqrt zwdth))
  (for ((= id 0) (< id ndim) (inc id))
    (def bdiff double
      (- (aref (fref (mref st-clt) bmax) id) (aref (fref (mref st-clt) bmin) id)))
    (if (< bdiff (* zeps zwdth))
        (begin
          (= (aref (fref (mref st-clt) bmax) id)
             (+ (aref (fref (mref st-clt) bmax) id)
                (* 0.5 (- (* zeps zwdth) bdiff))))
          (= (aref (fref (mref st-clt) bmin) id)
             (- (aref (fref (mref st-clt) bmin) id)
                (* 0.5 (- (* zeps zwdth) bdiff)))))))
  (= zwdth
     (* (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))
        (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))))
  (for ((= id 1) (< id ndim) (inc id))
    (= zwdth
       (* (- (aref (fref (mref st-clt) bmax) id) (aref (fref (mref st-clt) bmin) id))
          (- (aref (fref (mref st-clt) bmax) id)
             (aref (fref (mref st-clt) bmin) id)))))
  (= (fref (mref st-clt) zwdth) (csym::sqrt zwdth))
)

(def (csym::set-bndbox-cog st-clt zgmid nofc)
    (fn void (ptr cluster) (ptr (array double 3)) int) (decl ic int) (decl l int)
    (for ((= ic 0) (< ic (fref (mref st-clt) nnson)) (inc ic))
      (if (== ic 0) (begin (= l 0))
          (begin
            (= l
               (+ l (fref (mref (aref (fref (mref st-clt) pc-sons) (- ic 1))) nsize)))))
      (csym::set-bndbox-cog (aref (fref (mref st-clt) pc-sons) ic) (ptr (aref zgmid l)) nofc))
    (csym::cal-bndbox-cog st-clt zgmid nofc))

(def (csym::create-cluster nmbr ndpth nstrt nsize ndim nson)
    (fn (ptr cluster) int int int int int int)
  (decl st-clt (ptr cluster))
  (= st-clt (cast (ptr cluster) (csym::malloc (sizeof cluster)))) (= nmbr (+ nmbr 1))
  (= (fref (mref st-clt) nstrt) nstrt) (= (fref (mref st-clt) nsize) nsize)
  (= (fref (mref st-clt) ndim) ndim) (= (fref (mref st-clt) nnson) nson)
  (= (fref (mref st-clt) nmbr) nmbr) (= (fref (mref st-clt) ndpth) ndpth)
  ;;(= (fref (mref st-clt) pc-sons)
  ;;   (cast (ptr (ptr cluster)) (csym::malloc (* nson (sizeof (ptr cluster))))))
  (return st-clt)
)

(def (csym::free-st-clt st-clt) (fn void (ptr cluster))
  (decl ic int)
  (def nnson int (fref (mref st-clt) nnson))
  (for ((= ic 0) (< ic nnson) (inc ic))
    (free-st-clt (aref (fref (mref st-clt) pc-sons) ic)))
  (csym::free (fref (mref st-clt) bmin))
  (csym::free (fref (mref st-clt) bmax))
  (csym::free (fref (mref st-clt) pc-sons))
)

(def (csym::checkClusterTree f st-clt) (fn void (ptr FILE) (ptr cluster))
  (def width double (fref (mref st-clt) zwdth))
  ;;(csym::fprintf f "%lf~%" width)
  (if (== (fref (mref st-clt) nnson) 2)
    (begin 
      (csym::fprintf f "%lf~%" width)
      (csym::checkClusterTree f (aref (fref (mref st-clt) pc-sons) 0))
      (csym::checkClusterTree f (aref (fref (mref st-clt) pc-sons) 1))
    )
    (begin
      ;;(csym::fprintf f "%lf~%" width)
    )
  )
)

(def (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt nd md nclst)
    (wfn (ptr cluster) (ptr (array double 3)) (ptr (array double 3))
        int int int int int int)
  (decl i int)
  (decl id int)
  (decl il int)
  (decl nson int)
  (def minsz double (aref param 21))
  (def zcoef double (aref param 31))
  (decl zlmin (array double ndim))
  (decl zlmax (array double ndim))
  (def st-clt (ptr cluster))

  (= ndpth (+ ndpth 1))
  (if (<= nd minsz)
      (begin 
        (= nson 0)
        (= st-clt (csym::create-cluster nclst ndpth nsrt nd ndim nson)))
    (begin
      (if (> nd PN)
          (begin
            (def tsk (struct minmax))
            (minmax 0 nd zgmid (ptr tsk) ndpth nsrt)
            (= (aref zlmin 0) (fref tsk minx))
            (= (aref zlmax 0) (fref tsk maxx))
            (= (aref zlmin 1) (fref tsk miny))
            (= (aref zlmax 1) (fref tsk maxy))
            (= (aref zlmin 2) (fref tsk minz))
            (= (aref zlmax 2) (fref tsk maxz)))
	(begin
	  (for ((= id 0) (< id ndim) (inc id))
	    (= (aref zlmin id) (aref (aref zgmid 0) id))
	    (= (aref zlmax id) (aref zlmin id))
	    (for ((= il 1) (< il nd) (inc il))
	      (def zg double (aref (aref zgmid il) id))
	      (if (< zg (aref zlmin id)) (begin (= (aref zlmin id) zg))
		(if (< (aref zlmax id) zg) (begin (= (aref zlmax id) zg))))))))
      (def zdiff double (- (aref zlmax 0) (aref zlmin 0)))
      (def ncut int 0)
      (for ((= id 0) (< id ndim) (inc id))
	(def zidiff double (- (aref zlmax id) (aref zlmin id)))
	(if (> zidiff (* zcoef zdiff))
	    (begin
	      (= zdiff zidiff)
	      (= ncut id))))

      (def zlmid double (* 0.5 (+ (aref zlmax ncut) (aref zlmin ncut))))
      (def nl int 0) (def nr int (- nd 1))
      (if (> nd PN)
	  (begin
	    (def gn int (+ (/ nd CHUNK_SIZE) 1))
	    (def lessNum (ptr int))
	    (def moreNum (ptr int))
	    (def lessStart (ptr int))
	    (def moreStart (ptr int))
	    (= lessNum (cast (ptr int) (csym::malloc (* gn (sizeof int)))))
	    (= moreNum (cast (ptr int) (csym::malloc (* gn (sizeof int)))))
	    (= lessStart (cast (ptr int) (csym::malloc (* gn (sizeof int)))))
	    (= moreStart (cast (ptr int) (csym::malloc (* (+ gn 500000) (sizeof int)))))
	    (countNum 0 gn lessNum moreNum zgmid nd ncut zlmid gn ndpth nsrt)
	    (= nl (ssum 0 gn ndpth lessNum))
	    (= (aref lessStart 0) 0)
	    (= (aref moreStart 0) nl)
	    (if (and (!= nl 0) (!= nl nd))
		(begin 
		  (def tl int 0)
		  (def tm int nl)
		  (for ((= id 0) (< id (- gn 1)) (inc id))
		    (+= tl (aref lessNum id))
		    (+= tm (aref moreNum id))
		    (= (aref lessStart (+ id 1)) tl)
		    (= (aref moreStart (+ id 1)) tm)
		    )
		  (setStart 0 gn ncut gn nd zlmid lessNum moreNum lessStart moreStart zgmid tempzgmid ndpth nsrt)
		  ))
        (free lessNum)
	      (free moreNum)
	      (free lessStart)
	      (free moreStart)      
	    )
	(begin 
	  (while (< nl nr)
	    (while (and (< nl nd)
			(<= (aref (aref zgmid nl) ncut) zlmid))
	      (= nl (+ nl 1)))
	    (while (and (>= nr 0)
			(> (aref (aref zgmid nr) ncut) zlmid))
	      (= nr (- nr 1)))
	    (if (< nl nr)
		(begin
		  (for ((= id 0) (< id ndim) (inc id))
		    (def nh double (aref (aref zgmid nl) id))
		    (= (aref (aref zgmid nl) id) (aref (aref zgmid nr) id))
		    (= (aref (aref zgmid nr) id) nh)
		    )
		  )
	      )
	    )))
      (if (or (== nl nd) (== nl 0))
	  (begin
	    #|(= nson 1)
	    (= st-clt (csym::create-cluster nclst ndpth nsrt nd ndim nson))
	    (= (aref (fref (mref st-clt) pc-sons) 0)
	       (create-ctree-ssgeom (aref (fref (mref st-clt) pc-sons) 0) zgmid tempzgmid
					  ndpth ndscd nsrt nd md nclst)))|#
	    (= nson 0)
	    (= st-clt (csym::create-cluster nclst ndpth nsrt nd ndim nson)))
	(begin
	  (= nson 2)
	  (= st-clt (create-cluster nclst ndpth nsrt nd ndim nson))
    (= (fref (mref st-clt) pc-sons) (cast (ptr (ptr cluster)) (csym::malloc (* nson (sizeof (ptr cluster))))))
	  (def nsrt1 int nsrt)
	  (def nd1 int nl)
	  (def nsrt2 int (+ nsrt nl))
	  (def nd2 int (- nd nl))
	  (if (> nd PL)
	      (begin ;parallel here
		(if (<= nd PN)
		    (begin
		      (do-two
			  (= (aref (fref (mref st-clt) pc-sons) 0)
			     (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
			  (= (aref (fref (mref st-clt) pc-sons) 1)
			     (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
			(handles build-cluster
        ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
			  (:put
			    ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
			    (= this.zgmid (ptr (aref zgmid nl)))
			    (= this.tempzgmid (ptr (aref tempzgmid nl)))
			    (= this.ndpth ndpth)
			    (= this.ndscd ndscd)
			    (= this.nsrt nsrt2)
			    (= this.nd nd2)
			    (= this.md md)
			    (= this.nclst nclst))
			  (:get
			    (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
		      )
		  (begin
		    (if (bit-and ndpth 1)
			(begin ;look into next level
			  (if (and (> nd1 PN) (> nd2 PN))
			      (begin
				(do-two
				    (= (aref (fref (mref st-clt) pc-sons) 0)
				       (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
				    (= (aref (fref (mref st-clt) pc-sons) 1)
				       (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
				  (handles build-cluster
            ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
				    (:put
				      ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
				      (= this.zgmid (ptr (aref tempzgmid nl)))
				      (= this.tempzgmid (ptr (aref zgmid nl)))
				      (= this.ndpth ndpth)
				      (= this.ndscd ndscd)
				      (= this.nsrt nsrt2)
				      (= this.nd nd2)
				      (= this.md md)
				      (= this.nclst nclst)
				      )
				    (:get
				      (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
				)
			    (begin
			      (if (and (<= nd1 PN) (<= nd2 PN))
				  (begin
				    (for ((= i 0) (< i nd) (inc i))
				      (= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
				      (= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
				      (= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
				      )
				    (do-two
					(= (aref (fref (mref st-clt) pc-sons) 0)
					   (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
					(= (aref (fref (mref st-clt) pc-sons) 1)
					   (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
				      (handles build-cluster
              ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
					(:put
					  ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
					  (= this.zgmid (ptr (aref zgmid nl)))
					  (= this.tempzgmid (ptr (aref tempzgmid nl)))
					  (= this.ndpth ndpth)
					  (= this.ndscd ndscd)
					  (= this.nsrt nsrt2)
					  (= this.nd nd2)
					  (= this.md md)
					  (= this.nclst nclst))
					(:get
					  (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
				    )
				(begin
				  (if (and (<= nd1 PN) (> nd2 PN))
				      (begin
					(for ((= i 0) (< i nl) (inc i))
					  (= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
					  (= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
					  (= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
					  )
					(do-two
					    (= (aref (fref (mref st-clt) pc-sons) 0)
					       (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
					    (= (aref (fref (mref st-clt) pc-sons) 1)
					       (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
					  (handles build-cluster
              ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
					    (:put
					      ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
					      (= this.zgmid (ptr (aref tempzgmid nl)))
					      (= this.tempzgmid (ptr (aref zgmid nl)))
					      (= this.ndpth ndpth)
					      (= this.ndscd ndscd)
					      (= this.nsrt nsrt2)
					      (= this.nd nd2)
					      (= this.md md)
					      (= this.nclst nclst))
					    (:get
					      (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
					)

				    (begin
				      (if (and (> nd1 PN) (<= nd2 PN))
					  (begin
					    (for ((= i nl) (< i nd) (inc i))
					      (= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
					      (= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
					      (= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
					      )
					    (do-two
						(= (aref (fref (mref st-clt) pc-sons) 0)
						   (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
						(= (aref (fref (mref st-clt) pc-sons) 1)
						   (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
					      (handles build-cluster
                ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
						(:put
						  ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
						  (= this.zgmid (ptr (aref zgmid nl)))
						  (= this.tempzgmid (ptr (aref tempzgmid nl)))
						  (= this.ndpth ndpth)
						  (= this.ndscd ndscd)
						  (= this.nsrt nsrt2)
						  (= this.nd nd2)
						  (= this.md md)
						  (= this.nclst nclst))
						(:get
						  (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
					    )
					)
				      )
				    )
				  )
				
				)
			      )
			    )
			  )
		      (begin ;is even and exchange place
			(do-two
			    (= (aref (fref (mref st-clt) pc-sons) 0)
			       (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
			    (= (aref (fref (mref st-clt) pc-sons) 1)
			       (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
			  (handles build-cluster
        ;;(:guard (def p double 0)(if (> nd MB)(begin(= p 1))(begin(= p 0)))p)
			    (:put
			      ;;(= this.st-clt (aref (fref (mref st-clt) pc-sons) 1))
			      (= this.zgmid (ptr (aref tempzgmid nl)))
			      (= this.tempzgmid (ptr (aref zgmid nl)))
			      (= this.ndpth ndpth)
			      (= this.ndscd ndscd)
			      (= this.nsrt nsrt2)
			      (= this.nd nd2)
			      (= this.md md)
			      (= this.nclst nclst))
			    (:get
			      (= (aref (fref (mref st-clt) pc-sons) 1) this.r))))
			)
		      )
		    )
		  )
		)
	    (begin ;no parallel here
	      (if (<= nd PN)
		  (begin
		    (= (aref (fref (mref st-clt) pc-sons) 0)
		       (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
		    (= (aref (fref (mref st-clt) pc-sons) 1)
		       (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
		    )
		(begin
		  (if (bit-and ndpth 1)
		      (begin ;look into next level
			(if (and (> nd1 PN) (> nd2 PN))
			    (begin
			      (= (aref (fref (mref st-clt) pc-sons) 0)
				 (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
			      (= (aref (fref (mref st-clt) pc-sons) 1)
				 (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
			      )
			  (begin
			    (if (and (<= nd1 PN) (<= nd2 PN))
				(begin
				  (for ((= i 0) (< i nd) (inc i))
				    (= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
				    (= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
				    (= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
				    )
				  (= (aref (fref (mref st-clt) pc-sons) 0)
				     (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
				  (= (aref (fref (mref st-clt) pc-sons) 1)
				     (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
				  )
			      (begin
				(if (and (<= nd1 PN) (> nd2 PN))
				    (begin
				      (for ((= i 0) (< i nl) (inc i))
					(= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
					(= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
					(= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
					)
				      (= (aref (fref (mref st-clt) pc-sons) 0)
					 (create-ctree-ssgeom zgmid tempzgmid ndpth ndscd nsrt1 nd1 md nclst))
				      (= (aref (fref (mref st-clt) pc-sons) 1)
					 (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
				      )
				  (begin
				    (if (and (> nd1 PN) (<= nd2 PN))
					(begin
					  (for ((= i nl) (< i nd) (inc i))
					    (= (aref (aref zgmid i) 0) (aref (aref tempzgmid i) 0))
					    (= (aref (aref zgmid i) 1) (aref (aref tempzgmid i) 1))
					    (= (aref (aref zgmid i) 2) (aref (aref tempzgmid i) 2))
					    )
					  (= (aref (fref (mref st-clt) pc-sons) 0)
					     (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
					  (= (aref (fref (mref st-clt) pc-sons) 1)
					     (create-ctree-ssgeom (ptr (aref zgmid nl)) (ptr (aref tempzgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
					  )
				      )
				    )
				  )
				)
			      )
			    )
                          )
			)
		    (begin ;is even and exchange place
		      (= (aref (fref (mref st-clt) pc-sons) 0)
			 (create-ctree-ssgeom tempzgmid zgmid ndpth ndscd nsrt1 nd1 md nclst))
		      (= (aref (fref (mref st-clt) pc-sons) 1)
			 (create-ctree-ssgeom (ptr (aref tempzgmid nl)) (ptr (aref zgmid nl)) ndpth ndscd nsrt2 nd2 md nclst))
		      )
                    )
		  )
		)
	      )
	    )
	  )
	)))
  (= (fref (mref st-clt) ndscd) nd)
  ;;bounding box
  (= (fref (mref st-clt) bmin) (cast (ptr double) (csym::malloc (* 3 (sizeof double)))))
  (= (fref (mref st-clt) bmax) (cast (ptr double) (csym::malloc (* 3 (sizeof double)))))
  (def zeps double 1.0e-5)
  (if (> (fref (mref st-clt) nnson) 0)
    (begin
      (for ((= id 0) (< id ndim) (inc id))
        (= (aref (fref (mref st-clt) bmin) id) (aref (fref (mref (aref (fref (mref st-clt) pc-sons) 0)) bmin) id))
        (= (aref (fref (mref st-clt) bmax) id) (aref (fref (mref (aref (fref (mref st-clt) pc-sons) 0)) bmax) id))
      )
      (for ((= il 1) (< il (fref (mref st-clt) nnson)) (inc il))
        (for ((= id 0) (< id ndim) (inc id))
          (if (< (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmin) id) (aref (fref (mref st-clt) bmin) id))
            (begin
              (= (aref (fref (mref st-clt) bmin) id) (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmin) id))
            )
          )
          (if (< (aref (fref (mref st-clt) bmax) id) (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmax) id))
            (begin
              (= (aref (fref (mref st-clt) bmax) id) (aref (fref (mref (aref (fref (mref st-clt) pc-sons) il)) bmax) id))
            )
          )
        )
      )
    )
    (begin
      (for ((= id 0) (< id ndim) (inc id))
        (= (aref (fref (mref st-clt) bmin) id) (aref (aref zgmid 0) id))
        (= (aref (fref (mref st-clt) bmax) id) (aref (aref zgmid 0) id))
      )
      (for ((= id 0) (< id ndim) (inc id))
        (for ((= il 1) (< il (fref (mref st-clt) nsize)) (inc il))
          (if (< (aref (aref zgmid il) id) (aref (fref (mref st-clt) bmin) id))
            (begin
              (= (aref (fref (mref st-clt) bmin) id) (aref (aref zgmid il) id))
            )
          )
          (if (< (aref (fref (mref st-clt) bmax) id) (aref (aref zgmid il) id))
            (begin
              (= (aref (fref (mref st-clt) bmax) id) (aref (aref zgmid il) id))
            )
          )
        )
      )
    )
  )
  (def zwdth double (* (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0)) (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))))
  (for ((= id 1) (< id ndim) (inc id))
    (= zwdth
       (+ zwdth
          (* (- (aref (fref (mref st-clt) bmax) id)
                (aref (fref (mref st-clt) bmin) id))
             (- (aref (fref (mref st-clt) bmax) id)
                (aref (fref (mref st-clt) bmin) id)))))
  )
  (= zwdth (csym::sqrt zwdth))
  (for ((= id 0) (< id ndim) (inc id))
    (def bdiff double
      (- (aref (fref (mref st-clt) bmax) id) (aref (fref (mref st-clt) bmin) id)))
    (if (< bdiff (* zeps zwdth))
        (begin
          (= (aref (fref (mref st-clt) bmax) id)
             (+ (aref (fref (mref st-clt) bmax) id)
                (* 0.5 (- (* zeps zwdth) bdiff))))
          (= (aref (fref (mref st-clt) bmin) id)
             (- (aref (fref (mref st-clt) bmin) id)
                (* 0.5 (- (* zeps zwdth) bdiff)))))))
  (= zwdth
     (* (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))
        (- (aref (fref (mref st-clt) bmax) 0) (aref (fref (mref st-clt) bmin) 0))))
  (for ((= id 1) (< id ndim) (inc id))
    (= zwdth
       (* (- (aref (fref (mref st-clt) bmax) id) (aref (fref (mref st-clt) bmin) id))
          (- (aref (fref (mref st-clt) bmax) id)
             (aref (fref (mref st-clt) bmin) id)))))
  (= (fref (mref st-clt) zwdth) (csym::sqrt zwdth))
  ;;end of bounding box
  (return st-clt))

(def (minmax i1 i2 zgmid tsk ndpth nsrt) (wfn void int int (ptr (array double 3)) (ptr (struct minmax)) int int)
  ;;(csym::fprintf stderr "minmax: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID (+ i1 nsrt) (- (+ i2 nsrt) 1) ndpth)
  (def il int)
  (def minx double)
  (def miny double)
  (def minz double)
  (def maxx double)
  (def maxy double)
  (def maxz double)

  (= minx (aref (aref zgmid i1) 0))
  (= miny (aref (aref zgmid i1) 1))
  (= minz (aref (aref zgmid i1) 2))
  (= maxx minx)
  (= maxy miny)
  (= maxz minz)
  (do-many for il from (+ i1 1) to i2
    (if (== 0 (% (- i1 il 1) Check-interval))
	TCELL-CHECK-REQ)
    (def zg double)
    (= zg (aref (aref zgmid il) 0))
    (if (< zg minx) (begin (= minx zg))
      (if (< maxx zg) (begin (= maxx zg))))
    (= zg (aref (aref zgmid il) 1))
    (if (< zg miny) (begin (= miny zg))
      (if (< maxy zg) (begin (= maxy zg))))
    (= zg (aref (aref zgmid il) 2))
    (if (< zg minz) (begin (= minz zg))
      (if (< maxz zg) (begin (= maxz zg))))
    (handles minmax
      (:guard (< (- i2 il) 10000))
      (:put from j1 to j2
            (= this.i1 j1)
            (= this.i2 j2)
            (= this.zgmid zgmid)
            (= this.ndpth ndpth)
            (= this.nsrt nsrt))
      (:get
        (if (< this.minx minx) (begin (= minx this.minx)))
        (if (< maxx this.maxx) (begin (= maxx this.maxx)))
        (if (< this.miny miny) (begin (= miny this.miny)))
        (if (< maxy this.maxy) (begin (= maxy this.maxy)))
        (if (< this.minz minz) (begin (= minz this.minz)))
        (if (< maxz this.maxz) (begin (= maxz this.maxz))))))
  (= (fref tsk -> minx) minx)
  (= (fref tsk -> maxx) maxx)
  (= (fref tsk -> miny) miny)
  (= (fref tsk -> maxy) maxy)
  (= (fref tsk -> minz) minz)
  (= (fref tsk -> maxz) maxz)
)


(def (countNum i1 i2 lessNum moreNum zgmid nd ncut zlmid gn ndpth nsrt)
    (wfn void int int (ptr int) (ptr int) (ptr (array double 3)) int int double int int int)
  (def start int)
  (def end int)
  (= start (+ nsrt (* i1 CHUNK_SIZE)))
  (= end (- (+ nsrt (* i2 CHUNK_SIZE)) 1))
  (if (> end (+ nsrt nd))
    (begin (= end (- (+ nsrt nd) 1))
    )
  )
  (if (> start end) (begin (csym::fprintf stderr "error: start>end~%")))
  ;;(csym::fprintf stderr "countNum: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID start end ndpth)
  (def id int)
  (do-many for id from i1 to i2
    (if (== 0 (% (- i1 id 1) Check-interval))
	TCELL-CHECK-REQ)
    (def start int (* id CHUNK_SIZE))
    (def end int)
    (def im int)
    (def ln int 0)
    (def mn int 0)
    (if (== id (- gn 1))
	(begin (= end nd))
      (begin (= end (+ start CHUNK_SIZE))))
    (for ((= im start) (< im end) (inc im))
      (if (<= (aref (aref zgmid im) ncut) zlmid)
	  (begin (inc ln)) (begin (inc mn))))
    (= (aref lessNum id) ln)
    (= (aref moreNum id) mn)
    (handles countNum
      (:put from j1 to j2
	    (= this.i1 j1)
	    (= this.i2 j2)
	    (= this.lessNum lessNum)
	    (= this.moreNum moreNum)
	    (= this.ncut ncut)
	    (= this.zgmid zgmid)
	    (= this.nd nd)
	    (= this.gn gn)
	    (= this.zlmid zlmid)
	    (= this.ndpth ndpth)
      (= this.nsrt nsrt)
      )
      (:get)
      )
    )
  )

(def (setStart i1 i2 ncut gn nd zlmid lessNum moreNum lessStart moreStart zgmid tempzgmid ndpth nsrt)
  (wfn void int int int int int double (ptr int) (ptr int) (ptr int) (ptr int) (ptr (array double 3)) (ptr (array double 3)) int int)
  (def start int)
  (def end int)
  (= start (+ nsrt (* i1 CHUNK_SIZE)))
  (= end (- (+ nsrt (* i2 CHUNK_SIZE)) 1))
  (if (> end (+ nsrt nd))
    (begin (= end (- (+ nsrt nd) 1))
    )
  )
  (if (> start end) (begin (csym::fprintf stderr "error: start>end~%")))
  ;;(csym::fprintf stderr "setStart: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID start end ndpth)
  (def id int)
  (do-many for id from i1 to i2
    (if (== 0 (% (- i1 id 1) Check-interval))
        TCELL-CHECK-REQ)
    (def start int (* id CHUNK_SIZE))
    (def end int)
    (if (== id (- gn 1))
	(begin (= end nd))
      (begin (= end (+ start CHUNK_SIZE))))
    (def im int)
    (def ls int (aref lessStart id))
    (def ms int (aref moreStart id))
    (for ((= im start) (< im end) (inc im))
      (if (<= (aref (aref zgmid im) ncut) zlmid)
	  (begin
	    (= (aref (aref tempzgmid ls) 0) (aref (aref zgmid im) 0))
	    (= (aref (aref tempzgmid ls) 1) (aref (aref zgmid im) 1))
	    (= (aref (aref tempzgmid ls) 2) (aref (aref zgmid im) 2))
	    (inc ls)
	    )
	(begin
	  (= (aref (aref tempzgmid ms) 0) (aref (aref zgmid im) 0))
	  (= (aref (aref tempzgmid ms) 1) (aref (aref zgmid im) 1))
	  (= (aref (aref tempzgmid ms) 2) (aref (aref zgmid im) 2))
	  (inc ms))))
    (handles setStart
      (:put from j1 to j2
	    (= this.i1 j1)
	    (= this.i2 j2)
	    (= this.ncut ncut)
	    (= this.gn gn)
	    (= this.nd nd)
      (= this.nsrt nsrt)
	    (= this.zgmid zgmid)
	    (= this.zlmid zlmid)
	    (= this.tempzgmid tempzgmid)
	    (= this.lessNum lessNum)
	    (= this.moreNum moreNum)
	    (= this.lessStart lessStart)
	    (= this.moreStart moreStart)
      (= this.ndpth ndpth)
	)
      (:get)
	)
    )
  )

(def (ssum i1 i2 ndpth lessNum) (wfn int int int int (ptr int))
  ;;(csym::fprintf stderr "ssum: %f %d %d %d %d~%" (csym::get-wall-time) WORKER-ID i1 i2 ndpth)
  (def id int)
  (def sum int 0)
  (do-many for id from i1 to i2
    (if (== 0 (% (- i1 id 1) Check-interval))
        TCELL-CHECK-REQ)
    (+= sum (aref lessNum id))
    (handles ssum
      (:put from j1 to j2
	(= this.i1 j1)
	(= this.i2 j2)
	(= this.lessNum lessNum)
  (= this.ndpth ndpth)
	)
      (:get
	(+= sum this.sum)
	)
      )
    )
  (return sum)
)

(def (csym::get-wall-time) (fn double void)
  (decl time (struct timeval))
  (if (csym::gettimeofday (ptr time) NULL) (begin (return 0)))
  (return (+ (cast double (fref time tv-sec)) (* (cast double (fref time tv-usec)) 0.000001))))

(def (csym::get-cpu-time) (fn double void)
  (return (/ (cast double (csym::clock)) CLOCKS-PER-SEC)))