How to form a block diagonal matrix in Mata/ How to put blockdiag() in a loop?

[Joro Kolev]

Good afternoon,

I would like to form a block diagonal matrix which has square blocks of 1s on the diagonal, and the dimensions of these blocks are given in a Stata-vector. E.g.,

. sysuse auto, clear
(1978 automobile data)

. levelsof rep, matcell(Dimensions)
1 2 3 4 5

. mat list Dimensions

Dimensions[5,1]
    c1
r1   2
r2   8
r3  30
r4  18
r5  11

. return list

scalars:
                  r(N) =  69
                  r(r) =  5

macros:
             r(levels) : "1 2 3 4 5"

Now I want to form in Mata using blockdiag() a square block diagonal matrix of dimension 69 X 69, which has starting from the north-west corner of the matrix and going down to the south-east corner blocks of 1s on the main diagonal of dimenstions 2 X 2, 8 X 8, 30 X 30, 18 X 18, and 11 X 11.

If it is of any help, this is how I do this in Stata:

. mat Blocks = J(r(N),r(N),0)

. scalar J = 1

. forvalues l = 1/`r(r)' {
  2. mat Blocks[J,J] = J(Dimensions[`l',1], Dimensions[`l',1], 1)
  3. scalar J = J + Dimensions[`l',1]
  4. }

[John Mullahy]

Joro, might this work for your purposes

cap mata mata drop bdg()

set matastrict off

mata

function bdg(v) {
  real ma
  ma=J(v[1],v[1],1)
  for (i=2; i<=length(v); i++) {
    ma=blockdiag(ma,J(v[i],v[i],1))
  }
  return(ma)
}


end

Example

. mata bdg((2,3,5))
[symmetric]
         1    2    3    4    5    6    7    8    9   10
     +---------------------------------------------------+
   1 |   1                                               |
   2 |   1    1                                          |
   3 |   0    0    1                                     |
   4 |   0    0    1    1                                |
   5 |   0    0    1    1    1                           |
   6 |   0    0    0    0    0    1                      |
   7 |   0    0    0    0    0    1    1                 |
   8 |   0    0    0    0    0    1    1    1            |
   9 |   0    0    0    0    0    1    1    1    1       |
  10 |   0    0    0    0    0    1    1    1    1    1  |
     +---------------------------------------------------+

(added in edit)

Maybe this is a little more elegant?

function bdg(v) {
  real ma
  ma=J(0,0,.)
  for (i=1; i<=length(v); i++) {
    ma=blockdiag(ma,J(v[i],v[i],1))
  }
  return(ma)
}