/************************************************************************************ nwspgr: Nodes and weights for numerical integration on sparse grids (Smolyak) (c) Florian Heiss & Viktor Winschel Nov 05, 2007 Source code for compiling MATA functions *************************************************************************************/ clear version 9.2 mata: mata set matastrict on /************************************************************************************ nwspgr(): main function for generating nodes & weights for sparse grids intergration Syntax: nw = nwspgr(string scalar type, real scalar dim, real scalar k [, real scalar sym]) Input: type : type of 1D integration rule "KPU": Nested rule for unweighted integral over [0,1] "KPN": Nested rule for integral with Gaussian weight "GQU": Gaussian quadrature for unweighted integral over [0,1] (Gauss-Legendre) "GQN": Gaussian quadrature for integral with Gaussian weight (Gauss-Hermite) &func(): pointer to function func. func must accept level l and return matrix with the first column representing nodes and the second column representing weights. Declaration should look like: real matrix func(real scalar l) dim : dimension of the integration problem k : Accuracy level. The rule will be exact for polynomial up to total order 2k-1 sym : (optional) if sym==1 (or built-in 1D functions are used), the code is much faster with symmetric and incorrect with asymmetric 1D rules Output: nw : matrix with (dim+1) columns where the first dim columns represent the nodes and the last column represents the weights *************************************************************************************/ real matrix nwspgr(transmorphic scalar type, real scalar dim, real scalar k, | real scalar sym) { pointer colvector n1d, w1d real matrix nodes, newnw, is real colvector weights, sortvec, keep, R1d, Rq real rowvector midx real scalar minq, maxq, q, bq, j, lastkeep, r, nr, m, numnew, fh string scalar fn /* 1. get 1D nodes & weights (for positive orthant if symmetric) */ n1d = J(k,1,NULL) w1d = J(k,1,NULL) R1d = J(k,1,.) if (eltype(type)=="pointer") { for(q=1;q<=k;q++){ newnw = (*type)(q) numnew = rows(newnw) if (sym==1) newnw = newnw[floor(numnew/2)+1..numnew,.] R1d[q] = rows(newnw) n1d[q] = &(newnw[.,1]) w1d[q] = &(newnw[.,2]) } } else if (type=="GQU" | type=="GQN" | type=="KPU" | type=="KPN") { sym = 1 fn = findfile("nw"+type) fh = fopen(fn,"r") for(q=1;q<=k;q++){ newnw = fgetmatrix(fh) R1d[q] = rows(newnw) n1d[q] = &(newnw[.,1]) w1d[q] = &(newnw[.,2]) } fclose(fh) } else { errprintf("Unknown 1D integration rule") exit() } /* 1(b) initialization */ minq = max((0,k-dim)) maxq = k-1 nodes = J(0,dim,.) weights = J(0,1,.) /* 1(c) outer loop over q */ for(q=minq;q<=maxq;q++){ r = rows(weights) bq = (-1)^(maxq-q) * comb(dim-1,dim+q-k) /* matrix of all rowvectors in N^D_{D+q} */ is = SpGrGetSeq(dim,dim+q) Rq = R1d[is[.,1]] for(j=2;j<=dim;j++) Rq = Rq :* R1d[is[.,j]] nodes = (nodes \ J(colsum(Rq),dim,.)) weights = (weights \ J(colsum(Rq),1,.)) /* inner loop collecting product rules */ for (j=1;j<=rows(is);j++){ midx = is[j,.] newnw = SpGrKronProd(n1d[midx],w1d[midx]) nodes[r+1..r+Rq[j],.] = newnw[.,1..dim] weights[r+1..r+Rq[j],.] = bq:*newnw[.,dim+1] r=r+Rq[j] } /* collect equal nodes: first sort */ sortvec = order(nodes,1..dim) _collate(nodes,sortvec) _collate(weights,sortvec) /* then make a list of rows to keep and sum weights of equal nodes */ keep = 1; lastkeep = 1 for (r=2; r<=rows(nodes); r++){ if (nodes[r,.]==nodes[lastkeep,.]) weights[lastkeep]=weights[lastkeep]+weights[r] else{ lastkeep = r keep = (keep \ r) } } nodes = nodes[keep,.] weights = weights[keep] } // 2. expand rules to other orthants if symmetric rule if (sym==1){ nr = rows(nodes) m=*n1d[1] for(j=1;j<=dim;j++){ keep = J(nr,1,.) numnew = 0 for(r=1;r<=nr;r++) if (nodes[r,j] ~= m) { numnew++ keep[numnew]=r } if(numnew>0){ nodes = (nodes \ nodes[keep[1..numnew],.]) nodes[nr+1..nr+numnew,j] = 2*m :- nodes[nr+1..nr+numnew,j] weights = (weights \ weights[keep[1..numnew]]) nr=nr+numnew } } // 3. sort sortvec = order(nodes,1..dim) _collate(nodes,sortvec) _collate(weights,sortvec) } /* 4. return matrix of nodes & weights (normalization to account for rounding errors)*/ return((nodes, weights:/colsum(weights))) } /************************************************************************************ SpGrGetSeq(): function for generating matrix of all rowvectors in N^D_{norm} Syntax: out = nwspgr(d,norm) Input: d = dimension, will be #columns in output norm = row sum of elements each of the rows has to have Output: out = matrix with d columns. Each row represents one vector with all elements >=1 and the sum of elements == norm *************************************************************************************/ real matrix SpGrGetSeq(real scalar d, real scalar norm) { real rowvector seq real scalar a, fs, c seq = J(1,d,0) a=norm-d seq[1]=a fs = seq c=1 while (seq[d]