← Back to branch summary

~ubuntu-branches/ubuntu/raring/orpie/raring

« back to all changes in this revision

Viewing changes to gsl/mlgsl_linalg.c

  • Committer: Bazaar Package Importer
  • Author(s): Uwe Steinmann
  • Date: 2004-09-20 14:18:45 UTC
  • Revision ID: james.westby@ubuntu.com-20040920141845-j092sbrg4hd0nfsf
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
gsl/mlgsl_linalg.c
 
1
/* ocamlgsl - OCaml interface to GSL                        */
 
2
/* Copyright (�) 2002 - Olivier Andrieu                     */
 
3
/* distributed under the terms of the GPL version 2         */
 
4
 
 
5
 
 
6
#include <gsl/gsl_linalg.h>
 
7
 
 
8
#include "mlgsl_matrix_double.h"
 
9
#include "mlgsl_vector_double.h"
 
10
#include "mlgsl_permut.h"
 
11
 
 
12
 
 
13
/* simple matrix operations */
 
14
value ml_gsl_linalg_matmult_mod(value A, value omodA,
 
15
                                value B, value omodB,
 
16
                                value C)
 
17
{
 
18
  gsl_linalg_matrix_mod_t modA = Opt_arg(omodA, Int_val, GSL_LINALG_MOD_NONE);
 
19
  gsl_linalg_matrix_mod_t modB = Opt_arg(omodB, Int_val, GSL_LINALG_MOD_NONE);
 
20
  _DECLARE_MATRIX3(A, B, C);
 
21
  _CONVERT_MATRIX3(A, B, C);
 
22
  gsl_linalg_matmult_mod(&m_A, modA, &m_B, modB, &m_C);
 
23
  return Val_unit;
 
24
}
 
25
 
 
26
 
 
27
 
 
28
/* LU decomposition */
 
29
 
 
30
value ml_gsl_linalg_LU_decomp(value A, value P)
 
31
{
 
32
  int sign;
 
33
  GSL_PERMUT_OF_BIGARRAY(P);
 
34
  _DECLARE_MATRIX(A);
 
35
  _CONVERT_MATRIX(A);
 
36
  gsl_linalg_LU_decomp(&m_A, &perm_P, &sign);
 
37
  return Val_int(sign);
 
38
}
 
39
 
 
40
value ml_gsl_linalg_LU_solve(value LU, value P, value B, value X)
 
41
{
 
42
  GSL_PERMUT_OF_BIGARRAY(P);
 
43
  _DECLARE_MATRIX(LU);
 
44
  _DECLARE_VECTOR2(B,X);
 
45
  _CONVERT_MATRIX(LU);
 
46
  _CONVERT_VECTOR2(B,X);
 
47
  gsl_linalg_LU_solve(&m_LU, &perm_P, &v_B, &v_X);
 
48
  return Val_unit;
 
49
}
 
50
 
 
51
value ml_gsl_linalg_LU_svx(value LU, value P, value X)
 
52
{
 
53
  GSL_PERMUT_OF_BIGARRAY(P);
 
54
  _DECLARE_MATRIX(LU);
 
55
  _DECLARE_VECTOR(X);
 
56
  _CONVERT_MATRIX(LU);
 
57
  _CONVERT_VECTOR(X);
 
58
  gsl_linalg_LU_svx(&m_LU, &perm_P, &v_X);
 
59
  return Val_unit;
 
60
}
 
61
 
 
62
value ml_gsl_linalg_LU_refine(value A, value LU, value P, 
 
63
                              value B, value X, value RES)
 
64
{
 
65
  GSL_PERMUT_OF_BIGARRAY(P);
 
66
  _DECLARE_MATRIX2(A, LU);
 
67
  _DECLARE_VECTOR3(B, X, RES);
 
68
  _CONVERT_MATRIX2(A, LU);
 
69
  _CONVERT_VECTOR3(B, X, RES);
 
70
  gsl_linalg_LU_refine(&m_A, &m_LU, &perm_P, &v_B, &v_X, &v_RES);
 
71
  return Val_unit;
 
72
}
 
73
 
 
74
value ml_gsl_linalg_LU_refine_bc(value *argv, int argc)
 
75
{
 
76
  return ml_gsl_linalg_LU_refine(argv[0], argv[1], argv[2],
 
77
                                 argv[3], argv[4], argv[5]);
 
78
}
 
79
 
 
80
value ml_gsl_linalg_LU_invert(value LU, value P, value INV)
 
81
{
 
82
  GSL_PERMUT_OF_BIGARRAY(P);
 
83
  _DECLARE_MATRIX2(LU, INV);
 
84
  _CONVERT_MATRIX2(LU, INV);
 
85
  gsl_linalg_LU_invert(&m_LU, &perm_P, &m_INV);
 
86
  return Val_unit;
 
87
}
 
88
 
 
89
value ml_gsl_linalg_LU_det(value LU, value sig)
 
90
{
 
91
  _DECLARE_MATRIX(LU);
 
92
  _CONVERT_MATRIX(LU);
 
93
  return copy_double(gsl_linalg_LU_det(&m_LU, Int_val(sig)));
 
94
}
 
95
 
 
96
value ml_gsl_linalg_LU_lndet(value LU)
 
97
{
 
98
  _DECLARE_MATRIX(LU);
 
99
  _CONVERT_MATRIX(LU);
 
100
  return copy_double(gsl_linalg_LU_lndet(&m_LU));
 
101
}
 
102
 
 
103
value ml_gsl_linalg_LU_sgndet(value LU, value sig)
 
104
{
 
105
  _DECLARE_MATRIX(LU);
 
106
  _CONVERT_MATRIX(LU);
 
107
  return Val_int(gsl_linalg_LU_sgndet(&m_LU, Int_val(sig)));
 
108
}
 
109
 
 
110
 
 
111
 
 
112
/* QR decomposition */
 
113
value ml_gsl_linalg_QR_decomp(value A, value TAU)
 
114
{
 
115
  _DECLARE_MATRIX(A);
 
116
  _DECLARE_VECTOR(TAU);
 
117
  _CONVERT_MATRIX(A);
 
118
  _CONVERT_VECTOR(TAU);
 
119
  gsl_linalg_QR_decomp(&m_A, &v_TAU);
 
120
  return Val_unit;
 
121
}
 
122
 
 
123
 
 
124
value ml_gsl_linalg_QR_solve(value QR, value TAU, value B, value X)
 
125
{
 
126
  _DECLARE_MATRIX(QR);
 
127
  _DECLARE_VECTOR3(B,X,TAU);
 
128
  _CONVERT_MATRIX(QR);
 
129
  _CONVERT_VECTOR3(B,X,TAU);
 
130
  gsl_linalg_QR_solve(&m_QR, &v_TAU, &v_B, &v_X);
 
131
  return Val_unit;
 
132
}
 
133
 
 
134
value ml_gsl_linalg_QR_svx(value QR, value TAU, value X)
 
135
{
 
136
  _DECLARE_MATRIX(QR);
 
137
  _DECLARE_VECTOR2(TAU, X);
 
138
  _CONVERT_MATRIX(QR);
 
139
  _CONVERT_VECTOR2(TAU, X);
 
140
  gsl_linalg_QR_svx(&m_QR, &v_TAU, &v_X);
 
141
  return Val_unit;
 
142
}
 
143
 
 
144
value ml_gsl_linalg_QR_lssolve(value QR, value TAU, value B, value X, 
 
145
                               value RES)
 
146
{
 
147
  _DECLARE_MATRIX(QR);
 
148
  _DECLARE_VECTOR4(TAU, RES, B, X);
 
149
  _CONVERT_MATRIX(QR);
 
150
  _CONVERT_VECTOR4(TAU, RES, B, X);
 
151
  gsl_linalg_QR_lssolve(&m_QR, &v_TAU, &v_B, &v_X, &v_RES);
 
152
  return Val_unit;
 
153
}
 
154
 
 
155
value ml_gsl_linalg_QR_QTvec(value QR, value TAU, value V)
 
156
{
 
157
  _DECLARE_MATRIX(QR);
 
158
  _DECLARE_VECTOR2(TAU, V);
 
159
  _CONVERT_MATRIX(QR);
 
160
  _CONVERT_VECTOR2(TAU, V);
 
161
  gsl_linalg_QR_QTvec(&m_QR, &v_TAU, &v_V);
 
162
  return Val_unit;
 
163
}
 
164
 
 
165
value ml_gsl_linalg_QR_Qvec(value QR, value TAU, value V)
 
166
{
 
167
  _DECLARE_MATRIX(QR);
 
168
  _DECLARE_VECTOR2(TAU, V);
 
169
  _CONVERT_MATRIX(QR);
 
170
  _CONVERT_VECTOR2(TAU, V);
 
171
  gsl_linalg_QR_Qvec(&m_QR, &v_TAU, &v_V);
 
172
  return Val_unit;
 
173
}
 
174
 
 
175
value ml_gsl_linalg_QR_Rsolve(value QR, value B, value X)
 
176
{
 
177
  _DECLARE_MATRIX(QR);
 
178
  _DECLARE_VECTOR2(B,X);
 
179
  _CONVERT_MATRIX(QR);
 
180
  _CONVERT_VECTOR2(B,X);
 
181
  gsl_linalg_QR_Rsolve(&m_QR, &v_B, &v_X);
 
182
  return Val_unit;
 
183
}
 
184
 
 
185
value ml_gsl_linalg_QR_Rsvx(value QR, value X)
 
186
{
 
187
  _DECLARE_MATRIX(QR);
 
188
  _DECLARE_VECTOR(X);
 
189
  _CONVERT_MATRIX(QR);
 
190
  _CONVERT_VECTOR(X);
 
191
  gsl_linalg_QR_Rsvx(&m_QR, &v_X);
 
192
  return Val_unit;
 
193
}
 
194
 
 
195
value ml_gsl_linalg_QR_unpack(value QR, value TAU, value Q, value R)
 
196
{
 
197
  _DECLARE_MATRIX3(QR, Q, R);
 
198
  _DECLARE_VECTOR(TAU);
 
199
  _CONVERT_MATRIX3(QR, Q, R);
 
200
  _CONVERT_VECTOR(TAU);
 
201
  gsl_linalg_QR_unpack(&m_QR, &v_TAU, &m_Q, &m_R);
 
202
  return Val_unit;
 
203
}
 
204
 
 
205
value ml_gsl_linalg_QR_QRsolve(value Q, value R, value B, value X)
 
206
{
 
207
  _DECLARE_MATRIX2(Q, R);
 
208
  _DECLARE_VECTOR2(B, X);
 
209
  _CONVERT_MATRIX2(Q, R);
 
210
  _CONVERT_VECTOR2(B, X);
 
211
  gsl_linalg_QR_QRsolve(&m_Q, &m_R, &v_B, &v_X);
 
212
  return Val_unit;
 
213
}
 
214
 
 
215
value ml_gsl_linalg_QR_update(value Q, value R, value W, value V)
 
216
{
 
217
  _DECLARE_MATRIX2(Q, R);
 
218
  _DECLARE_VECTOR2(W, V);
 
219
  _CONVERT_MATRIX2(Q, R);
 
220
  _CONVERT_VECTOR2(W, V);
 
221
  gsl_linalg_QR_update(&m_Q, &m_R, &v_W, &v_V);
 
222
  return Val_unit;
 
223
}
 
224
 
 
225
value ml_gsl_linalg_R_solve(value R, value B, value X)
 
226
{
 
227
  _DECLARE_MATRIX(R);
 
228
  _DECLARE_VECTOR2(B, X);
 
229
  _CONVERT_MATRIX(R);
 
230
  _CONVERT_VECTOR2(B, X);
 
231
  gsl_linalg_R_solve(&m_R, &v_B, &v_X);
 
232
  return Val_unit;
 
233
}
 
234
 
 
235
/* missing ? */
 
236
/*  value ml_gsl_linalg_R_svx(value R, value X) */
 
237
/*  { */
 
238
/*    DECLARE_MATRIX(R); */
 
239
/*    DECLARE_VECTOR(X); */
 
240
/*    gsl_linalg_R_svx(&m_R, &v_X); */
 
241
/*    return Val_unit; */
 
242
/*  } */
 
243
 
 
244
 
 
245
 
 
246
/* QR Decomposition with Column Pivoting */
 
247
value ml_gsl_linalg_QRPT_decomp(value A, value TAU, value P, value NORM)
 
248
{
 
249
  int signum;
 
250
  GSL_PERMUT_OF_BIGARRAY(P);
 
251
  _DECLARE_MATRIX(A);
 
252
  _DECLARE_VECTOR2(TAU, NORM);
 
253
  _CONVERT_MATRIX(A);
 
254
  _CONVERT_VECTOR2(TAU, NORM);
 
255
  gsl_linalg_QRPT_decomp(&m_A, &v_TAU, &perm_P, &signum, &v_NORM);
 
256
  return Val_int(signum);
 
257
}
 
258
 
 
259
value ml_gsl_linalg_QRPT_decomp2(value A, value Q, value R,
 
260
                                 value TAU, value P, value NORM)
 
261
{
 
262
  int signum;
 
263
  GSL_PERMUT_OF_BIGARRAY(P);
 
264
  _DECLARE_MATRIX3(A, Q, R);
 
265
  _DECLARE_VECTOR2(TAU, NORM);
 
266
  _CONVERT_MATRIX3(A, Q, R);
 
267
  _CONVERT_VECTOR2(TAU, NORM);
 
268
  gsl_linalg_QRPT_decomp2(&m_A, &m_Q, &m_R, &v_TAU, &perm_P, &signum, &v_NORM);
 
269
  return Val_int(signum);
 
270
}
 
271
 
 
272
value ml_gsl_linalg_QRPT_decomp2_bc(value *argv, int argc)
 
273
{
 
274
  return ml_gsl_linalg_QRPT_decomp2(argv[0], argv[1], argv[2],
 
275
                                    argv[3], argv[4], argv[5]);
 
276
}
 
277
 
 
278
value ml_gsl_linalg_QRPT_solve(value QR, value TAU, value P, 
 
279
                               value B, value X)
 
280
{
 
281
  GSL_PERMUT_OF_BIGARRAY(P);
 
282
  _DECLARE_MATRIX(QR);
 
283
  _DECLARE_VECTOR3(TAU, B, X);
 
284
  _CONVERT_MATRIX(QR);
 
285
  _CONVERT_VECTOR3(TAU, B, X);
 
286
  gsl_linalg_QRPT_solve(&m_QR, &v_TAU, &perm_P, &v_B, &v_X);
 
287
  return Val_unit;
 
288
}
 
289
 
 
290
value ml_gsl_linalg_QRPT_svx(value QR, value TAU, value P, value X)
 
291
{
 
292
  GSL_PERMUT_OF_BIGARRAY(P);
 
293
  _DECLARE_MATRIX(QR);
 
294
  _DECLARE_VECTOR2(TAU, X);
 
295
  _CONVERT_MATRIX(QR);
 
296
  _CONVERT_VECTOR2(TAU, X);
 
297
  gsl_linalg_QRPT_svx(&m_QR, &v_TAU, &perm_P, &v_X);
 
298
  return Val_unit;
 
299
}
 
300
 
 
301
value ml_gsl_linalg_QRPT_QRsolve(value Q, value R, value P, 
 
302
                                 value B, value X)
 
303
{
 
304
  GSL_PERMUT_OF_BIGARRAY(P);
 
305
  _DECLARE_MATRIX2(Q, R);
 
306
  _DECLARE_VECTOR2(B, X);
 
307
  _CONVERT_MATRIX2(Q, R);
 
308
  _CONVERT_VECTOR2(B, X);
 
309
  gsl_linalg_QRPT_QRsolve(&m_Q, &m_R, &perm_P, &v_B, &v_X);
 
310
  return Val_unit;
 
311
}
 
312
 
 
313
value ml_gsl_linalg_QRPT_update(value Q, value R, value P, 
 
314
                                value U, value V)
 
315
{
 
316
  GSL_PERMUT_OF_BIGARRAY(P);
 
317
  _DECLARE_MATRIX2(Q, R);
 
318
  _DECLARE_VECTOR2(U, V);
 
319
  _CONVERT_MATRIX2(Q, R);
 
320
  _CONVERT_VECTOR2(U, V);
 
321
  gsl_linalg_QRPT_update(&m_Q, &m_R, &perm_P, &v_U, &v_V);
 
322
  return Val_unit;
 
323
}
 
324
 
 
325
value ml_gsl_linalg_QRPT_Rsolve(value QR, value P, 
 
326
                                value B, value X)
 
327
{
 
328
  GSL_PERMUT_OF_BIGARRAY(P);
 
329
  _DECLARE_MATRIX(QR);
 
330
  _DECLARE_VECTOR2(B, X);
 
331
  _CONVERT_MATRIX(QR);
 
332
  _CONVERT_VECTOR2(B, X);
 
333
  gsl_linalg_QRPT_Rsolve(&m_QR, &perm_P, &v_B, &v_X);
 
334
  return Val_unit;
 
335
}
 
336
 
 
337
value ml_gsl_linalg_QRPT_Rsvx(value QR, value P, value X)
 
338
{
 
339
  GSL_PERMUT_OF_BIGARRAY(P);
 
340
  _DECLARE_MATRIX(QR);
 
341
  _DECLARE_VECTOR(X);
 
342
  _CONVERT_MATRIX(QR);
 
343
  _CONVERT_VECTOR(X);
 
344
  gsl_linalg_QRPT_Rsvx(&m_QR, &perm_P, &v_X);
 
345
  return Val_unit;
 
346
}
 
347
 
 
348
 
 
349
 
 
350
/* Singular Value Decomposition */
 
351
value ml_gsl_linalg_SV_decomp(value A, value V, value S, value WORK)
 
352
{
 
353
  _DECLARE_MATRIX2(A, V);
 
354
  _DECLARE_VECTOR2(S, WORK);
 
355
  _CONVERT_MATRIX2(A, V);
 
356
  _CONVERT_VECTOR2(S, WORK);
 
357
  gsl_linalg_SV_decomp(&m_A, &m_V, &v_S, &v_WORK);
 
358
  return Val_unit;
 
359
}
 
360
 
 
361
value ml_gsl_linalg_SV_decomp_mod(value A, value X, value V, 
 
362
                                  value S, value WORK)
 
363
{
 
364
  _DECLARE_MATRIX3(A, V, X);
 
365
  _DECLARE_VECTOR2(S, WORK);
 
366
  _CONVERT_MATRIX3(A, V, X);
 
367
  _CONVERT_VECTOR2(S, WORK);
 
368
  gsl_linalg_SV_decomp_mod(&m_A, &m_X, &m_V, &v_S, &v_WORK);
 
369
  return Val_unit;
 
370
}
 
371
 
 
372
value ml_gsl_linalg_SV_decomp_jacobi(value A, value V, value S)
 
373
{
 
374
  _DECLARE_MATRIX2(A, V);
 
375
  _DECLARE_VECTOR(S);
 
376
  _CONVERT_MATRIX2(A, V);
 
377
  _CONVERT_VECTOR(S);
 
378
  gsl_linalg_SV_decomp_jacobi(&m_A, &m_V, &v_S);
 
379
  return Val_unit;
 
380
}
 
381
 
 
382
value ml_gsl_linalg_SV_solve(value U, value V, value S, value B, value X)
 
383
{
 
384
  _DECLARE_MATRIX2(U, V);
 
385
  _DECLARE_VECTOR3(S, B, X);
 
386
  _CONVERT_MATRIX2(U, V);
 
387
  _CONVERT_VECTOR3(S, B, X);
 
388
  gsl_linalg_SV_solve(&m_U, &m_V, &v_S, &v_B, &v_X);
 
389
  return Val_unit;
 
390
}
 
391
 
 
392
 
 
393
 
 
394
/* Cholesky decomposition */
 
395
 
 
396
value ml_gsl_linalg_cholesky_decomp(value A)
 
397
{
 
398
  _DECLARE_MATRIX(A);
 
399
  _CONVERT_MATRIX(A);
 
400
  gsl_linalg_cholesky_decomp(&m_A);
 
401
  return Val_unit;
 
402
}
 
403
 
 
404
value ml_gsl_linalg_cholesky_solve(value CHO, value B, value X)
 
405
{
 
406
  _DECLARE_MATRIX(CHO);
 
407
  _DECLARE_VECTOR2(B, X);
 
408
  _CONVERT_MATRIX(CHO);
 
409
  _CONVERT_VECTOR2(B, X);
 
410
  gsl_linalg_cholesky_solve(&m_CHO, &v_B, &v_X);
 
411
  return Val_unit;
 
412
}
 
413
 
 
414
value ml_gsl_linalg_cholesky_svx(value CHO, value X)
 
415
{
 
416
  _DECLARE_MATRIX(CHO);
 
417
  _DECLARE_VECTOR(X);
 
418
  _CONVERT_MATRIX(CHO);
 
419
  _CONVERT_VECTOR(X);
 
420
  gsl_linalg_cholesky_svx(&m_CHO, &v_X);
 
421
  return Val_unit;
 
422
}
 
423
 
 
424
 
 
425
 
 
426
/* Tridiagonal Decomposition of Real Symmetric Matrices */
 
427
value ml_gsl_linalg_symmtd_decomp(value A, value TAU)
 
428
{
 
429
  _DECLARE_MATRIX(A);
 
430
  _DECLARE_VECTOR(TAU);
 
431
  _CONVERT_MATRIX(A);
 
432
  _CONVERT_VECTOR(TAU);
 
433
  gsl_linalg_symmtd_decomp(&m_A, &v_TAU);
 
434
  return Val_unit;
 
435
}
 
436
 
 
437
value ml_gsl_linalg_symmtd_unpack(value A, value TAU, value Q, 
 
438
                                  value DIAG, value SUBDIAG)
 
439
{
 
440
  _DECLARE_MATRIX2(A, Q);
 
441
  _DECLARE_VECTOR3(TAU, DIAG, SUBDIAG);
 
442
  _CONVERT_MATRIX2(A, Q);
 
443
  _CONVERT_VECTOR3(TAU, DIAG, SUBDIAG);
 
444
  gsl_linalg_symmtd_unpack(&m_A, &v_TAU, &m_Q, &v_DIAG, &v_SUBDIAG);
 
445
  return Val_unit;
 
446
}
 
447
 
 
448
value ml_gsl_linalg_symmtd_unpack_T(value A, value DIAG, value SUBDIAG)
 
449
{
 
450
  _DECLARE_MATRIX(A);
 
451
  _DECLARE_VECTOR2(DIAG, SUBDIAG);
 
452
  _CONVERT_MATRIX(A);
 
453
  _CONVERT_VECTOR2(DIAG, SUBDIAG);
 
454
  gsl_linalg_symmtd_unpack_T(&m_A, &v_DIAG, &v_SUBDIAG);
 
455
  return Val_unit;
 
456
}
 
457
 
 
458
 
 
459
 
 
460
/* Tridiagonal Decomposition of Hermitian Matrices */
 
461
 
 
462
 
 
463
 
 
464
/* Bidiagonalization */
 
465
value ml_gsl_linalg_bidiag_decomp(value A, value TAU_U, value TAU_V)
 
466
{
 
467
  _DECLARE_MATRIX(A);
 
468
  _DECLARE_VECTOR2(TAU_U, TAU_V);
 
469
  _CONVERT_MATRIX(A);
 
470
  _CONVERT_VECTOR2(TAU_U, TAU_V);
 
471
  gsl_linalg_bidiag_decomp(&m_A, &v_TAU_U, &v_TAU_V);
 
472
  return Val_unit;
 
473
}
 
474
 
 
475
value ml_gsl_linalg_bidiag_unpack(value A, value TAU_U, value U, 
 
476
                                  value TAU_V, value V, 
 
477
                                  value DIAG, value SUPERDIAG)
 
478
{
 
479
  _DECLARE_MATRIX3(A, U, V);
 
480
  _DECLARE_VECTOR4(TAU_U, TAU_V, DIAG, SUPERDIAG);
 
481
  _CONVERT_MATRIX3(A, U, V);
 
482
  _CONVERT_VECTOR4(TAU_U, TAU_V, DIAG, SUPERDIAG);
 
483
  gsl_linalg_bidiag_unpack(&m_A, &v_TAU_U, &m_U, &v_TAU_V, &m_V,
 
484
                           &v_DIAG, &v_SUPERDIAG);
 
485
  return Val_unit;
 
486
}
 
487
 
 
488
value ml_gsl_linalg_bidiag_unpack_bc(value *argv, int argc)
 
489
{
 
490
  return ml_gsl_linalg_bidiag_unpack(argv[0], argv[1], argv[2], 
 
491
                                     argv[3], argv[4], argv[5], argv[6]);
 
492
}
 
493
 
 
494
value ml_gsl_linalg_bidiag_unpack2(value A, value TAU_U, value TAU_V, value V)
 
495
{
 
496
  _DECLARE_MATRIX2(A, V);
 
497
  _DECLARE_VECTOR2(TAU_U, TAU_V);
 
498
  _CONVERT_MATRIX2(A, V);
 
499
  _CONVERT_VECTOR2(TAU_U, TAU_V);
 
500
  gsl_linalg_bidiag_unpack2(&m_A, &v_TAU_U, &v_TAU_V, &m_V);
 
501
  return Val_unit;
 
502
}
 
503
 
 
504
value ml_gsl_linalg_bidiag_unpack_B(value A, value DIAG, value SUPERDIAG)
 
505
{
 
506
  _DECLARE_MATRIX(A);
 
507
  _DECLARE_VECTOR2(DIAG, SUPERDIAG);
 
508
  _CONVERT_MATRIX(A);
 
509
  _CONVERT_VECTOR2(DIAG, SUPERDIAG);
 
510
  gsl_linalg_bidiag_unpack_B(&m_A, &v_DIAG, &v_SUPERDIAG);
 
511
  return Val_unit;
 
512
}
 
513
 
 
514
 
 
515
 
 
516
/* Householder solver */
 
517
value ml_gsl_linalg_HH_solve(value A, value B, value X)
 
518
{
 
519
  _DECLARE_MATRIX(A);
 
520
  _DECLARE_VECTOR2(B,X);
 
521
  _CONVERT_MATRIX(A);
 
522
  _CONVERT_VECTOR2(B,X);
 
523
  gsl_linalg_HH_solve(&m_A, &v_B, &v_X);
 
524
  return Val_unit;
 
525
}
 
526
 
 
527
value ml_gsl_linalg_HH_svx(value A, value X)
 
528
{
 
529
  _DECLARE_MATRIX(A);
 
530
  _DECLARE_VECTOR(X);
 
531
  _CONVERT_MATRIX(A);
 
532
  _CONVERT_VECTOR(X);
 
533
  gsl_linalg_HH_svx(&m_A, &v_X);
 
534
  return Val_unit;
 
535
}
 
536
 
 
537
 
 
538
/* Tridiagonal Systems */
 
539
value ml_gsl_linalg_solve_symm_tridiag(value DIAG, value E, value B, value X)
 
540
{
 
541
  _DECLARE_VECTOR4(DIAG, E, B, X);
 
542
  _CONVERT_VECTOR4(DIAG, E, B, X);
 
543
  gsl_linalg_solve_symm_tridiag(&v_DIAG, &v_E, &v_B, &v_X);
 
544
  return Val_unit;
 
545
}
 
546
 
 
547
value ml_gsl_linalg_solve_tridiag(value DIAG, value ABOVE, value BELOW,
 
548
                                  value B, value X)
 
549
{
 
550
  _DECLARE_VECTOR5(DIAG, ABOVE, BELOW, B, X);
 
551
  _CONVERT_VECTOR5(DIAG, ABOVE, BELOW, B, X);
 
552
  gsl_linalg_solve_tridiag(&v_DIAG, &v_ABOVE, &v_BELOW, &v_B, &v_X);
 
553
  return Val_unit;
 
554
}
 
555
 
 
556
value ml_gsl_linalg_solve_symm_cyc_tridiag(value DIAG, value E, value B, value X)
 
557
{
 
558
  _DECLARE_VECTOR4(DIAG, E, B, X);
 
559
  _CONVERT_VECTOR4(DIAG, E, B, X);
 
560
  gsl_linalg_solve_symm_cyc_tridiag(&v_DIAG, &v_E, &v_B, &v_X);
 
561
  return Val_unit;
 
562
}
 
563
 
 
564
value ml_gsl_linalg_solve_cyc_tridiag(value DIAG, value ABOVE, value BELOW,
 
565
                                      value B, value X)
 
566
{
 
567
  _DECLARE_VECTOR5(DIAG, ABOVE, BELOW, B, X);
 
568
  _CONVERT_VECTOR5(DIAG, ABOVE, BELOW, B, X);
 
569
  gsl_linalg_solve_cyc_tridiag(&v_DIAG, &v_ABOVE, &v_BELOW, &v_B, &v_X);
 
570
  return Val_unit;
 
571
}
 
572
 
 
573
 
 
574
/* exponential */
 
575
#define GSL_MODE_val Int_val
 
576
 
 
577
value ml_gsl_linalg_exponential_ss(value A, value eA, value mode)
 
578
{
 
579
  _DECLARE_MATRIX2(A, eA);
 
580
  _CONVERT_MATRIX2(A, eA);
 
581
  gsl_linalg_exponential_ss(&m_A, &m_eA, GSL_MODE_val(mode));
 
582
  return Val_unit;
 
583
}