@@ -9,41 +9,64 @@ static char *imcinname;
9
9
static char * outimname ;
10
10
static char * outcoeffname ;
11
11
12
- static CLICMDARGDEF farg [] = {
13
- {CLIARG_IMG , ".inc" , "input 3D cube" , "imc" , CLIARG_VISIBLE_DEFAULT , (void * * )& imcinname , NULL },
14
- {CLIARG_STR , ".outm" , "output modes" , "outm" , CLIARG_VISIBLE_DEFAULT , (void * * )& outimname , NULL },
15
- {CLIARG_STR , ".outcoeff" , "output coeffs" , "outcoeff" , CLIARG_VISIBLE_DEFAULT , (void * * )& outcoeffname , NULL }};
16
-
17
- static CLICMDDATA CLIcmddata = {"imsvd" , "Singular values decomposition" , CLICMD_FIELDS_DEFAULTS };
12
+ static CLICMDARGDEF farg [] = {{CLIARG_IMG ,
13
+ ".inc" ,
14
+ "input 3D cube" ,
15
+ "imc" ,
16
+ CLIARG_VISIBLE_DEFAULT ,
17
+ (void * * ) & imcinname ,
18
+ NULL },
19
+ {CLIARG_STR ,
20
+ ".outm" ,
21
+ "output modes" ,
22
+ "outm" ,
23
+ CLIARG_VISIBLE_DEFAULT ,
24
+ (void * * ) & outimname ,
25
+ NULL },
26
+ {CLIARG_STR ,
27
+ ".outcoeff" ,
28
+ "output coeffs" ,
29
+ "outcoeff" ,
30
+ CLIARG_VISIBLE_DEFAULT ,
31
+ (void * * ) & outcoeffname ,
32
+ NULL }};
33
+
34
+ static CLICMDDATA CLIcmddata = {
35
+ "imsvd" , "Singular values decomposition" , CLICMD_FIELDS_DEFAULTS };
18
36
19
37
// detailed help
20
- static errno_t help_function () { return RETURN_SUCCESS ; }
38
+ static errno_t help_function ()
39
+ {
40
+ return RETURN_SUCCESS ;
41
+ }
21
42
22
43
// rotation matrix written as SVD_VTm
23
44
24
- errno_t linopt_compute_SVDdecomp (const char * IDin_name , const char * IDout_name , const char * IDcoeff_name ,
25
- imageID * outID )
45
+ errno_t linopt_compute_SVDdecomp (const char * IDin_name ,
46
+ const char * IDout_name ,
47
+ const char * IDcoeff_name ,
48
+ imageID * outID )
26
49
{
27
50
DEBUG_TRACE_FSTART ();
28
51
29
- imageID IDin ;
30
- imageID IDout ;
31
- imageID IDcoeff ;
32
- gsl_matrix * matrix_D ; /* input */
33
- gsl_matrix * matrix_Dtra ;
34
- gsl_matrix * matrix_DtraD ;
35
- gsl_matrix * matrix_DtraD_evec ;
36
- gsl_vector * matrix_DtraD_eval ;
52
+ imageID IDin ;
53
+ imageID IDout ;
54
+ imageID IDcoeff ;
55
+ gsl_matrix * matrix_D ; /* input */
56
+ gsl_matrix * matrix_Dtra ;
57
+ gsl_matrix * matrix_DtraD ;
58
+ gsl_matrix * matrix_DtraD_evec ;
59
+ gsl_vector * matrix_DtraD_eval ;
37
60
gsl_eigen_symmv_workspace * w ;
38
- gsl_matrix * matrix_save ;
61
+ gsl_matrix * matrix_save ;
39
62
40
- long m ;
41
- long n ;
63
+ long m ;
64
+ long n ;
42
65
uint32_t * arraysizetmp ;
43
66
44
67
imageID ID_VTmatrix ;
45
68
46
- arraysizetmp = (uint32_t * )malloc (sizeof (uint32_t ) * 3 );
69
+ arraysizetmp = (uint32_t * ) malloc (sizeof (uint32_t ) * 3 );
47
70
if (arraysizetmp == NULL )
48
71
{
49
72
FUNC_RETURN_FAILURE ("malloc returns NULL pointer" );
@@ -58,33 +81,44 @@ errno_t linopt_compute_SVDdecomp(const char *IDin_name, const char *IDout_name,
58
81
m = data .image [IDin ].md [0 ].size [2 ];
59
82
60
83
matrix_DtraD_eval = gsl_vector_alloc (m );
61
- matrix_D = gsl_matrix_alloc (n , m );
62
- matrix_Dtra = gsl_matrix_alloc (m , n );
63
- matrix_DtraD = gsl_matrix_alloc (m , m );
84
+ matrix_D = gsl_matrix_alloc (n , m );
85
+ matrix_Dtra = gsl_matrix_alloc (m , n );
86
+ matrix_DtraD = gsl_matrix_alloc (m , m );
64
87
matrix_DtraD_evec = gsl_matrix_alloc (m , m );
65
88
66
89
/* write matrix_D */
67
90
for (long k = 0 ; k < m ; k ++ )
68
91
{
69
92
for (long ii = 0 ; ii < n ; ii ++ )
70
93
{
71
- gsl_matrix_set (matrix_D , ii , k , data .image [IDin ].array .F [k * n + ii ]);
94
+ gsl_matrix_set (matrix_D ,
95
+ ii ,
96
+ k ,
97
+ data .image [IDin ].array .F [k * n + ii ]);
72
98
}
73
99
}
74
100
/* compute DtraD */
75
- gsl_blas_dgemm (CblasTrans , CblasNoTrans , 1.0 , matrix_D , matrix_D , 0.0 , matrix_DtraD );
101
+ gsl_blas_dgemm (CblasTrans ,
102
+ CblasNoTrans ,
103
+ 1.0 ,
104
+ matrix_D ,
105
+ matrix_D ,
106
+ 0.0 ,
107
+ matrix_DtraD );
76
108
77
109
/* compute the inverse of DtraD */
78
110
79
111
/* first, compute the eigenvalues and eigenvectors */
80
- w = gsl_eigen_symmv_alloc (m );
112
+ w = gsl_eigen_symmv_alloc (m );
81
113
matrix_save = gsl_matrix_alloc (m , m );
82
114
gsl_matrix_memcpy (matrix_save , matrix_DtraD );
83
115
gsl_eigen_symmv (matrix_save , matrix_DtraD_eval , matrix_DtraD_evec , w );
84
116
85
117
gsl_matrix_free (matrix_save );
86
118
gsl_eigen_symmv_free (w );
87
- gsl_eigen_symmv_sort (matrix_DtraD_eval , matrix_DtraD_evec , GSL_EIGEN_SORT_ABS_DESC );
119
+ gsl_eigen_symmv_sort (matrix_DtraD_eval ,
120
+ matrix_DtraD_evec ,
121
+ GSL_EIGEN_SORT_ABS_DESC );
88
122
89
123
create_2Dimage_ID (IDcoeff_name , m , 1 , & IDcoeff );
90
124
@@ -96,22 +130,33 @@ errno_t linopt_compute_SVDdecomp(const char *IDin_name, const char *IDout_name,
96
130
/** Write rotation matrix to go from DM modes to eigenmodes */
97
131
arraysizetmp [0 ] = m ;
98
132
arraysizetmp [1 ] = m ;
99
- ID_VTmatrix = image_ID ("SVD_VTm" );
133
+ ID_VTmatrix = image_ID ("SVD_VTm" );
100
134
if (ID_VTmatrix != -1 )
101
135
{
102
136
delete_image_ID ("SVD_VTm" , DELETE_IMAGE_ERRMODE_WARNING );
103
137
}
104
- create_image_ID ("SVD_VTm" , 2 , arraysizetmp , _DATATYPE_FLOAT , 0 , 0 , 0 , & ID_VTmatrix );
138
+ create_image_ID ("SVD_VTm" ,
139
+ 2 ,
140
+ arraysizetmp ,
141
+ _DATATYPE_FLOAT ,
142
+ 0 ,
143
+ 0 ,
144
+ 0 ,
145
+ & ID_VTmatrix );
105
146
for (long ii = 0 ; ii < m ; ii ++ ) // modes
106
147
for (long k = 0 ; k < m ; k ++ ) // modes
107
148
{
108
- data .image [ID_VTmatrix ].array .F [k * m + ii ] = (float )gsl_matrix_get (matrix_DtraD_evec , k , ii );
149
+ data .image [ID_VTmatrix ].array .F [k * m + ii ] =
150
+ (float ) gsl_matrix_get (matrix_DtraD_evec , k , ii );
109
151
}
110
152
111
153
/// Compute SVD decomp
112
154
113
- FUNC_CHECK_RETURN (create_3Dimage_ID (IDout_name , data .image [IDin ].md [0 ].size [0 ], data .image [IDin ].md [0 ].size [1 ],
114
- data .image [IDin ].md [0 ].size [2 ], & IDout ));
155
+ FUNC_CHECK_RETURN (create_3Dimage_ID (IDout_name ,
156
+ data .image [IDin ].md [0 ].size [0 ],
157
+ data .image [IDin ].md [0 ].size [1 ],
158
+ data .image [IDin ].md [0 ].size [2 ],
159
+ & IDout ));
115
160
116
161
for (long kk = 0 ; kk < m ; kk ++ ) /// eigen mode index
117
162
{
@@ -122,7 +167,8 @@ errno_t linopt_compute_SVDdecomp(const char *IDin_name, const char *IDout_name,
122
167
for (long ii = 0 ; ii < n ; ii ++ )
123
168
{
124
169
data .image [IDout ].array .F [kk * n + ii ] +=
125
- data .image [ID_VTmatrix ].array .F [kk1 * m + kk ] * data .image [IDin ].array .F [kk1 * n + ii ];
170
+ data .image [ID_VTmatrix ].array .F [kk1 * m + kk ] *
171
+ data .image [IDin ].array .F [kk1 * n + ii ];
126
172
}
127
173
}
128
174
}
0 commit comments