Restructuring a GSL matrix

I have a set of K*K square matrices and I would like to create a new (K+1)*(K+1) square matrices where each K size block matrix copies in this new one then there would be a row and column of zeros except its K+1-th element which would be a constant C. The final blocks would be all zeros as well. Here there is :

1
2
3
4
5
6
7
8
9
10
  array([[ 0.09990096, -0.03676744, -0.04634088, -0.03676744,  0.01353258,  0.01704259, -0.04634088,  0.01704259,  0.02123895],
       [-0.03676744,  0.07944992, -0.03132864,  0.01353258, -0.0292484 ,  0.01152807,  0.01704259, -0.03700603,  0.01447335],
       [-0.04634088, -0.03132864,  0.13579896,  0.01704259,  0.01152807, -0.04996734,  0.02123895,  0.01447335, -0.06275539],
       [-0.03676744,  0.01353258,  0.01704259,  0.07944992, -0.0292484 , -0.03700603, -0.03132864,  0.01152807,  0.01447335],
       [ 0.01353258, -0.0292484 ,  0.01152807, -0.0292484 ,  0.06309598, -0.02495025,  0.01152807, -0.02495025,  0.00981151],
       [ 0.01704259,  0.01152807, -0.04996734, -0.03700603, -0.02495025,  0.1081392 ,  0.01447335,  0.00981151, -0.04253148],
       [-0.04634088,  0.01704259,  0.02123895, -0.03132864,  0.01152807,  0.01447335,  0.13579896, -0.04996734, -0.06275539],
       [ 0.01704259, -0.03700603,  0.01447335,  0.01152807, -0.02495025,  0.00981151, -0.04996734,  0.1081392 , -0.04253148],
       [ 0.02123895,  0.01447335, -0.06275539,  0.01447335,  0.00981151, -0.04253148, -0.06275539, -0.04253148,  0.18437719]])

new matrix should look like this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
array([[ 0.09990096, -0.03676744, -0.04634088,  0.        , -0.03676744,   0.01353258,  0.01704259,  0.        , -0.04634088,  0.01704259,   0.02123895,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [-0.03676744,  0.07944992, -0.03132864,  0.        ,  0.01353258,   -0.0292484 ,  0.01152807, 0.        ,  0.01704259, -0.03700603,   0.01447335,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [-0.04634088, -0.03132864,  0.13579896,  0.        ,  0.01704259,   0.01152807, -0.04996734,  0.        ,  0.02123895,  0.01447335,  -0.06275539,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  5.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [-0.03676744,  0.01353258,  0.01704259,  0.        ,  0.07944992,  -0.0292484 , -0.03700603,  0.        , -0.03132864,  0.01152807,   0.01447335,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.01353258, -0.0292484 ,  0.01152807,  0.        , -0.0292484 ,   0.06309598, -0.02495025,  0.        ,  0.01152807, -0.02495025,   0.00981151,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.01704259,  0.01152807, -0.04996734,  0.        , -0.03700603,  -0.02495025,  0.1081392 ,  0.        ,  0.01447335,  0.00981151,  -0.04253148,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  5.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [-0.04634088,  0.01704259,  0.02123895,  0.        , -0.03132864,   0.01152807,  0.01447335,  0.        ,  0.13579896, -0.04996734,  -0.06275539,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.01704259, -0.03700603,  0.01447335,  0.        ,  0.01152807,  -0.02495025,  0.00981151,  0.        , -0.04996734,  0.1081392 ,  -0.04253148,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.02123895,  0.01447335, -0.06275539,  0.        ,  0.01447335,   0.00981151, -0.04253148,  0.        , -0.06275539, -0.04253148,   0.18437719,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  5.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  5.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  5.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  5.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,   0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  5.        ]])

How can I do it using GSL_matrix?
Your matrix is python (but missing the 'd'), not C++.
I'm assuming you want a C++ solution.

So, K in your example is 3? And you want to expand every 3x3 "block" into a 4x4 block, ultimately creating a 16x16 matrix, filling in the new entries with 0's except for the bottom-left corner of the block, which will be a constanct C (5 in your example). Any extra blocks needed to fill out the full 16x16 matrix are all 0's.

Seems simple enough. Where are you stuck?
I want the c++ answer! I wrote this code just the place of constant values are messed up right now.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <stdio.h>
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

int main () {
double Qx[]={ 0.09990096, -0.03676744, -0.04634088, -0.03676744,  0.01353258,  0.01704259, -0.04634088,  0.01704259,  0.02123895,
             -0.03676744,  0.07944992, -0.03132864,  0.01353258, -0.0292484 ,  0.01152807,  0.01704259, -0.03700603,  0.01447335,
             -0.04634088, -0.03132864,  0.13579896,  0.01704259,  0.01152807, -0.04996734,  0.02123895,  0.01447335, -0.06275539,
             -0.03676744,  0.01353258,  0.01704259,  0.07944992, -0.0292484 , -0.03700603, -0.03132864,  0.01152807,  0.01447335,
              0.01353258, -0.0292484 ,  0.01152807, -0.0292484 ,  0.06309598, -0.02495025,  0.01152807, -0.02495025,  0.00981151,
              0.01704259,  0.01152807, -0.04996734, -0.03700603, -0.02495025,  0.1081392 ,  0.01447335,  0.00981151, -0.04253148,
             -0.04634088,  0.01704259,  0.02123895, -0.03132864,  0.01152807,  0.01447335,  0.13579896, -0.04996734, -0.06275539,
              0.01704259, -0.03700603,  0.01447335,  0.01152807, -0.02495025,  0.00981151, -0.04996734,  0.1081392 , -0.04253148,
              0.02123895,  0.01447335, -0.06275539,  0.01447335,  0.00981151, -0.04253148, -0.06275539, -0.04253148,  0.18437719};
    int nr=3;

    gsl_matrix_view Qmatrix = gsl_matrix_view_array (Qx, nr*nr, nr*nr);
    gsl_matrix *Qexp= gsl_matrix_calloc((nr+1)*(nr+1),(nr+1)*(nr+1));
    gsl_matrix_view Qmatrix_view;
    gsl_matrix_view Qexp_view;
    double cons=5.;
    for (int row=0; row<nr; ++row){
        for (int col=0; col<nr; ++col){
        Qmatrix_view = gsl_matrix_submatrix(&Qmatrix.matrix, row*nr, col*nr,nr,nr);
        Qexp_view    = gsl_matrix_submatrix(Qexp,row*(nr+1),col*(nr+1),nr,nr);
        gsl_matrix_set(Qexp, (row+1)*nr,(col+1)*nr, cons);
        gsl_matrix_memcpy(&Qexp_view.matrix, &Qmatrix_view.matrix);
            
        }
    }
    printf("\n---\n---\n---\n");
    for (int row=0; row<(nr+1)*(nr+1); ++row)
        for (int col=0; col<(nr+1)*(nr+1); ++col)
            printf(col==nr*(nr+1)?"%6.5f\n":"%6.5f ",gsl_matrix_get(Qexp,row,col));
    printf("\n---\n---\n");
    return 0;
}
Last edited on
I'm not familiar with gsl, but taking a guess, try moving gsl_matrix_set to after gsl_matrix_memcpy.

EDIT: Note that the manual says that for gsl_matrix_memcpy:
The two matrices must have the same size.

EDIT2: Oh wait, they are the same size. Oh well. I guess I can't see the error.
Last edited on
Topic archived. No new replies allowed.