Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: mata: why this loop is so slow


From   "U.A.QU.AP" <uaquap@gnet.tn>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: mata: why this loop is so slow
Date   Wed, 31 May 2006 00:22:17 +0100 (CET)

I'm really in trouble my computer ( centrino 1.7 GHz 512 RAM)
run Bill's Code in 1.84 second.
But when changing slightly Bill's code (namely dimension of
matrix q and h) in order to match exactly my problem execution
time took about 16 MINUTES see log file above. Any comments
. do "C:\projets\nicolas\toulouse\gould.do"

. cscript
-------------------------------------------------------------------------BEGIN

r; t=0.00 23:35:02

. 
. set memory 30m
(30720k)
r; t=0.00 23:35:02

. set rmsg on
r; t=0.00 23:35:02

. local OBS 35900
r; t=0.00 23:35:02

. 
. set obs `OBS'
obs was 0, now 35900
r; t=0.00 23:35:02

. set seed 1234
r; t=0.00 23:35:02

. 
. forvalues i=1(1)46 {
  2.        gen v`i' = uniform()
  3. }
r; t=0.41 23:35:03

. forvalues i=1(1)35 {
  2.        gen vprobit`i' = uniform()
  3. }
r; t=0.20 23:35:03

. 
. gen dmills  = uniform()
r; t=0.00 23:35:03

. gen selection = uniform()
r; t=0.02 23:35:03

. gen t = mod(_n-1, 5) + 1
r; t=0.02 23:35:03

. 
. local instrume v1  v2  v3  v4  v5  v6  v7  v8  v9  v10 ///
>               v11 v12 v13 v14 v15 v16 v17 v18 v19 v20 ///
>               v21 v22 v23 v24 v25 v26 v27 v28 v29 v30 ///
>               v31 v32 v33 v34 v35 v36 v37 v38 v39 v40 ///
>               v41 v42 v43 v44 v45 v46
r; t=0.00 23:35:03

. 
. local vprobit vprobit1 vprobit2 vprobit3 vprobit4 vprobit5
vprobit6 vprobit7 vprobit8 vprobit9 vprobit10 
> ///
>                 vprobit11 vprobit12 vprobit13 vprobit14
vprobit15 vprobit16 vprobit17 vprobit18 vprobit19
>  vprobit20 ///
>                 vprobit21 vprobit22 vprobit23 vprobit24
vprobit25 vprobit26 vprobit27 vprobit28 vprobit29
>  vprobit30 ///
>                 vprobit31 vprobit32 vprobit33 vprobit34
vprobit35
r; t=0.00 23:35:03

. mata:
------------------------------------------------- mata (type
end to exit) ---------------------------------

: w = st_data(., tokens(st_local("instrume")))    // matrix

: q = st_data(., tokens(st_local("vprobit")))     // matrix

: dmills=st_data(., "dmills")                     // colvector

: s=st_data(., "selection")                       // colvector

: t=st_data(., "t")                               // colvector

: 
: chw=cols(w)                                     // scalar

: cq=cols(q)                                      // scalar

: 
: h = uniform(`OBS', 70)

: 
: n = nt = `OBS'

: T  = 21

: beta = uniform(46,1)

: 
: F=J(cols(h),(cols(q)*T),0)

: gg=cols(F)

: bb=rows(F)

: gg
  735

: bb
  70

: for (i=1; i<=nt; i++) {
>        jacobi=J(chw, (cq*T) ,0)
>        jacobi[(chw - T + t[i]), (cq*t[i] -cq+1)..cq*t[i]] =
>                                                  
q[i,.]:*dmills[i,.]
>        F = F + (s[i,.]:*h[i,.]')*beta'*jacobi
> }

: 
: F=(1/n)*F

: gg=cols(F)

: bb=rows(F)

: gg 
  735

: bb
  70

: 
: end
-----------------------------------------------------------------------------------------------------------
r; t=974.58 23:51:17

. 
end of do-file
r; t=975.23 23:51:17

. 



---- Original message ---- 
>Date: Tue, 30 May 2006 12:41:43 -0500
>From: wgould@stata.com (William Gould, Stata)  
>Subject: Re: st: mata:  why this loop is so slow  
>To: statalist@hsphsun2.harvard.edu
>
>AbdelRahmen <uaquap@gnet.tn> reported
>
>> I've written this portion of code (see below) and it works
>> perfectly the problem is the execution time which is about 15
>> mn. Under SAS (IML) my co-author code the same thing and the
>> execution time is 1 (one) minute on 35900 observations. is
>> there a problem with my programmation Style?
>
>I replied, AbdelRahmen replied, and I now have run his code
on my 
>computer.  
>
>It ran is 0.84 seconds on my computer.
>
>Perhaps I have done something wrong.  I concocted phony
values, but 
>everything I think is as AbdelRahmen describes it.  I wonder how 
>fast the attached do-file runs on his computer.
>
>-- Bill
>wgould@stata.com
>
>File try.do follows:
>------------------------------------------------------------------------------
>cscript
>
>set memory 30m
>set rmsg on
>local OBS 35900
>
>set obs `OBS'
>set seed 1234
>
>forvalues i=1(1)46 {
>        gen v`i' = uniform()
>}
>gen vprobit = uniform()
>gen dmills  = uniform()
>gen selection = uniform()
>gen t = mod(_n-1, 5) + 1
>
>local instrume v1  v2  v3  v4  v5  v6  v7  v8  v9  v10 ///
>               v11 v12 v13 v14 v15 v16 v17 v18 v19 v20 ///
>               v21 v22 v23 v24 v25 v26 v27 v28 v29 v30 ///
>               v31 v32 v33 v34 v35 v36 v37 v38 v39 v40 ///
>               v41 v42 v43 v44 v45 v46 
>
>local vprobit "vprobit"
>
>mata:
>
>w = st_data(., tokens(st_local("instrume")))    // matrix
>q = st_data(., tokens(st_local("vprobit")))     // matrix
>dmills=st_data(., "dmills")                     // colvector
>s=st_data(., "selection")                       // colvector
>t=st_data(., "t")                               // colvector
>
>chw=cols(w)                                     // scalar
>cq=cols(q)                                      // scalar
>
>h = uniform(`OBS', 1)
>
>n = nt = `OBS'
>T  = 21
>beta = uniform(46,1)
>
>F = J(1,21,0)
>for (i=1; i<=nt; i++) {
>        jacobi=J(chw, (cq*T) ,0)
>        jacobi[(chw - T + t[i]), (cq*t[i] -cq+1)..cq*t[i]] = 
>                                                  
q[i,.]:*dmills[i,.]
>        F = F + (s[i,.]:*h[i,.]')*beta'*jacobi
>}
>
>F=(1/n)*F
>
>end
>------------------------------------------------------------------------------
><end>
>*
>*   For searches and help try:
>*   http://www.stata.com/support/faqs/res/findit.html
>*   http://www.stata.com/support/statalist/faq
>*   http://www.ats.ucla.edu/stat/stata/
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index