/*
Filename: mylib_matrix.mac
Contents: private maxima library on matrix
2022/12/8: created by Hideo Kodama  
2022/12/15: version 1.00
2023/4/21: version 1.10
Last update: 2023/4/21
*/

/******>Top*/

/*****> Matrix functions */
/****>proc: matrixtype */
/* matrixtype(M) => [rownum, colnum]
*/
matrixtype(M) := block(
    [m,n],
    m: length(listme(col(M,1))),
    n: length(listme(row(M,1))),
    return ([m,n])
)$ /* end proc:matrixtpe*/

/****>proc: blockMatrices */
/* blockMatrices(mxlist) 
   mxlist = [M1, ..., Mn] => M : blockdiagonal
*/
blockMatrices(mxlist) := block( [ a,b, j0,k0, rownl, colnl, M],
    mxlist : map(lambda([x], if not matrixp(x) then ( matrix([x])) else (x) ), mxlist),
    rownl : [], colnl: [],
    for x in mxlist do (
        [a, b] : matrix_size (x),
        rownl : append( rownl, [a] ),
        colnl : append(colnl, [b])
    ), /*end do*/
    M : zeromatrix( sumlist(rownl), sumlist(colnl)),
    [j0, k0] : [0,0],
    for i: 1 thru length(mxlist) do (
        for j: 1 thru rownl[i] do (
            for k: 1 thru colnl[i] do (
                M[j0+j, k0+k] : mxlist[i][j,k]
            )  /*end do*/
        ), /*end do*/
        j0 : j0+rownl[i],  k0 : k0+colnl[i]
    ), /*end do*/
    return (copy(M))
)$  /*end proc: blockMatrices */

/****>proc: colMatrices */
/* colMatrices( mxlist) 
   mxlist = [M1, ..., Mn], M1 =matrix(a1, b), ..., Mn = matrix(an,b)
   => M = matrix(a1+...+an,b)
*/
colMatrices( mxlist ) := block( [x, n, rown, coln, M],
    n : length(mxlist), 
    x : matrix_size(mxlist[1]),
    coln : x[2], 
    M : mxlist[1],
    for i:2 thru n do (
       x : matrix_size(mxlist[i]),
       if x[2]# coln then (
          error("colMatrices: all the component matrices should have the same column number")
       ),
       rown : x[1],
       M : apply(addrow, cons(M, makelist( row(mxlist[i],j), j, 1, rown) ) )
   ), /*end do*/
   return (M)
)$ /*end proc: colMatrices */


/****>proc: rowMatrices */
/* rowMatrices( mxlist) 
   mxlist = [M1, ..., Mn], M1 =matrix(a, b1), ..., Mn = matrix(a,bn)
   => M = matrix(a,b1+...+bn)
*/
rowMatrices( mxlist ) := block( [x, n, rown, coln, M],
    n : length(mxlist), 
    x : matrix_size(mxlist[1]),
    rown : x[1], 
    M : mxlist[1],
    for i:2 thru n do (
       x : matrix_size(mxlist[i]),
       if x[1]# rown then (
          error("rowMatrices: all the component matrices should have the same number of rows")
       ),
       coln : x[2],
       M : apply(addcol, cons(M, makelist( col(mxlist[i],j), j, 1, coln) ) )
   ), /*end do*/
   return (M)
)$ /*end proc: rowMatrices */

/****>proc: submatrix1 */
/* submatrix1( rowrange, M, colrange ) => M1 : submatrix of M
   rowrange = [i_1, i_p],
   colrange = [j_1, j_q]
   M1 = Matrix(p,q) 
*/
submatrix1 (rowrange, M, colrange) := block ( [m, n, rowl, coll, p, q, M1],
    [m, n] : matrix_size(M),
    rowl: makelist( i, i, rowrange[1], rowrange[2]),
    coll: makelist( i, i, colrange[1], colrange[2]),
    p : length(rowl), q : length(coll),
    M1 : zeromatrix(p,q),
    for i : 1 thru p do (
        for j:1 thru q do (
            M1[i,j] : M[rowl[i], coll[j]]
        )  /*end do*/
    ), /*end do*/
    return (copy(M1))
)$   /*end proc: submatrix1 */

/****>proc:listme */
/* listme = list_matrix_entries */
listme ( M ) := block(
    return ( list_matrix_entries (M))
)$  /*end proc: listme */

/******>End */