Removing a column from a gsl_matrix

Question

I would like to remove one column of a gsl_matrix (the i-th column) and then copy its content into another gsl_matrix. This is my code:

#include <stdio.h>
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

int remove_col (int K,
                int N,
                int i,//between range 1 to N
                gsl_matrix *Sn, //Kx(N-1)
                gsl_matrix *Z)
{
    gsl_matrix_view Z_view;

    gsl_matrix_view Sn_view;
    if (i==0){
       Z_view  = gsl_matrix_submatrix (Z, 0, 1, K, N-1);
       gsl_matrix_memcpy (Sn, &Z_view.matrix);
    }else{
       Z_view  = gsl_matrix_submatrix (Z, 0, 0, K, i);
       Sn_view = gsl_matrix_submatrix (Sn, 0, 0, K, i);
       gsl_matrix_memcpy (&Sn_view.matrix, &Z_view.matrix);
       Z_view  = gsl_matrix_submatrix (Z, 0, i, K, N-i);
       Sn_view = gsl_matrix_submatrix (Sn, 0, i-1, K, N-i);
       gsl_matrix_memcpy (&Sn_view.matrix, &Z_view.matrix);
    }     
    for ( int row = 0; row < K; ++row ) { 
        for ( int col = 0; col < N-1; ++col ) { 
            printf( "\t%3.1f",  gsl_matrix_get( Sn , row, col));
        }
        printf( "\n" );
    };
}
int main () {
double Zx[]={1, 1, 1, 1, 1, 1,
             0, 1, 1, 1, 1, 1,
             0, 0, 1, 1, 1, 1,
             0, 1, 0, 1, 1, 1,
             1, 1, 0, 1, 1, 0,
             1, 1, 1, 1, 1, 1,
             1, 1, 0, 1, 0, 0,
             1, 1, 1, 0, 1, 0,
             1, 1, 0, 1, 0, 0,
             0, 1, 1, 1, 0, 1};

int r=10;             
int c=6; 
gsl_matrix_view Z = gsl_matrix_view_array (Zx, r, c); 

gsl_matrix *R= gsl_matrix_calloc(c,r);
for ( int row = 0; row < r; ++row ) { 
           for ( int col = 0; col < c; ++col ) { 
               gsl_matrix_set( R, col, row, gsl_matrix_get( &Z.matrix , row, col));
           }
}

gsl_matrix *Sn = gsl_matrix_alloc(c,r-1);
remove_col (c,r,1, Sn, R);
for (int row=0; row<c; ++row)      
    for (int col=0; col<r-1; ++col)  
        printf(col==r-2?"%6.3f\n":"%6.3f ",gsl_matrix_get(Sn,row,col));
gsl_matrix_free(R);    
gsl_matrix_free(Sn);    
} 
    return 0;
}

I get this error:

gsl: submatrix_source.c:37: ERROR: first dimension overflows matrix
Default GSL error handler invoked.
Aborted (core dumped)

Is my approach correct? Any suggestion?

Answer 1

This version works:

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

int remove_col (int K,
                int N,
                int i,//between range 1 to N
                gsl_matrix *Sn, //Kx(N-1)
                gsl_matrix *Z) {
  gsl_matrix_view Z_view;

  gsl_matrix_view Sn_view;
  if (i==0){
    Z_view  = gsl_matrix_submatrix (Z, 0, 1, K, N-1);
    gsl_matrix_memcpy (Sn, &Z_view.matrix);
  } else{
    Z_view  = gsl_matrix_submatrix (Z, 0, 0, K, i);
    Sn_view = gsl_matrix_submatrix (Sn, 0, 0, K, i);
    gsl_matrix_memcpy (&Sn_view.matrix, &Z_view.matrix);
    Z_view  = gsl_matrix_submatrix (Z, 0, i, K, N-i);
    Sn_view = gsl_matrix_submatrix (Sn, 0, i-1, K, N-i);
    gsl_matrix_memcpy (&Sn_view.matrix, &Z_view.matrix);
  }     
}  
int main () {
  double Zx[]={1, 1, 1, 1, 1, 1,
               0, 1, 1, 1, 1, 1,
               0, 0, 1, 1, 1, 1,
               0, 1, 0, 1, 1, 1,
               1, 1, 0, 1, 1, 0,
               1, 1, 1, 1, 1, 1,
               1, 1, 0, 1, 0, 0,
               1, 1, 1, 0, 1, 0,
               1, 1, 0, 1, 0, 0,
               0, 1, 1, 1, 0, 1};

  int r=10;             
  int c=6; 
  gsl_matrix_view Z = gsl_matrix_view_array (Zx, r, c); 

  gsl_matrix *R= gsl_matrix_calloc(c,r);
  for ( int row = 0; row < r; ++row ) { 
    for ( int col = 0; col < c; ++col ) { 
      gsl_matrix_set( R, col, row, gsl_matrix_get( &Z.matrix , row, col));
    }
  }
  for (int row=0; row<r; ++row){
    for (int col=0; col<c; ++col)
      printf(col==r-2?"%6.3f\n":"%6.3f ",gsl_matrix_get(R,col,row));
    printf("\n");
  }

  gsl_matrix *Sn = gsl_matrix_alloc(c,r-1);
  remove_col (c,r,1, Sn, R);
  for (int row=0; row<r-1; ++row){
    for (int col=0; col<c; ++col)
      printf(col==r-2?"%6.3f\n":"%6.3f ",gsl_matrix_get(Sn,col,row));
    printf("\n");
  }
  gsl_matrix_free(R);    
  gsl_matrix_free(Sn);    

  return 0;
}