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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include "ltfat.h"
#include "ltfat_types.h"
LTFAT_EXTERN void
LTFAT_NAME(gabtight_fac)(const LTFAT_COMPLEX *gf, const ltfatInt L,const ltfatInt R,
const ltfatInt a, const ltfatInt M,
LTFAT_COMPLEX *gtightf)
{
ltfatInt h_a, h_m;
LTFAT_COMPLEX *Sf, *U, *VT, *gfwork;
LTFAT_REAL *S;
const LTFAT_COMPLEX zzero = (LTFAT_COMPLEX) 0.0;//{0.0, 0.0 };
const LTFAT_COMPLEX alpha = (LTFAT_COMPLEX) (1.0+0.0*I);//{1.0, 0.0 };
const ltfatInt N=L/a;
const ltfatInt c=gcd(a, M,&h_a, &h_m);
const ltfatInt p=a/c;
const ltfatInt q=M/c;
const ltfatInt d=N/q;
S = ltfat_malloc(p*sizeof*S);
Sf = ltfat_malloc(p*p*sizeof*Sf);
U = ltfat_malloc(p*p*sizeof*U);
VT = ltfat_malloc(p*q*R*sizeof*VT);
gfwork = ltfat_malloc(L*R*sizeof*gfwork);
/* Copy the contents of gf to gfwork because LAPACK overwrites
* the input.
*/
memcpy(gfwork,gf,L*R*sizeof*gfwork);
for (ltfatInt rs=0; rs<c*d; rs++)
{
/* Compute the thin SVD */
LTFAT_NAME(ltfat_gesvd)(p, q*R, gfwork+rs*p*q*R, p,
S, U, p, VT, p);
/* Combine U and V. */
LTFAT_NAME(ltfat_gemm)(CblasNoTrans,CblasNoTrans,p,q*R,p,
&alpha,(const LTFAT_COMPLEX*)U,p,
(const LTFAT_COMPLEX*)VT,p,
&zzero,gtightf+rs*p*q*R, p);
}
LTFAT_SAFEFREEALL(gfwork,Sf,S,U,VT);
}
LTFAT_EXTERN void
LTFAT_NAME(gabtightreal_fac)(const LTFAT_COMPLEX *gf, const ltfatInt L, const ltfatInt R,
const ltfatInt a, const ltfatInt M,
LTFAT_COMPLEX *gtightf)
{
ltfatInt h_a, h_m;
LTFAT_COMPLEX *Sf, *U, *VT, *gfwork;
LTFAT_REAL *S;
const LTFAT_COMPLEX zzero = (LTFAT_COMPLEX) 0.0;
const LTFAT_COMPLEX alpha = (LTFAT_COMPLEX) 1.0+0.0*I;//{1.0, 0.0 };
const ltfatInt N=L/a;
const ltfatInt c=gcd(a, M,&h_a, &h_m);
const ltfatInt p=a/c;
const ltfatInt q=M/c;
const ltfatInt d=N/q;
/* This is a floor operation. */
const ltfatInt d2= d/2+1;
S = ltfat_malloc(p*sizeof*S);
Sf = ltfat_malloc(p*p*sizeof*Sf);
U = ltfat_malloc(p*p*sizeof*U);
VT = ltfat_malloc(p*q*R*sizeof*VT);
gfwork = ltfat_malloc(L*R*sizeof*gfwork);
/* Copy the contents of gf to gfwork because LAPACK overwrites
* the input.
*/
memcpy(gfwork,gf,L*R*sizeof*gfwork);
for (ltfatInt rs=0; rs<c*d2; rs++)
{
/* Compute the thin SVD */
LTFAT_NAME(ltfat_gesvd)(p, q*R, gfwork+rs*p*q*R, p,
S, U, p, VT, p);
/* Combine U and V. */
LTFAT_NAME(ltfat_gemm)(CblasNoTrans,CblasNoTrans,p,q*R,p,
&alpha,(const LTFAT_COMPLEX*)U,p,
(const LTFAT_COMPLEX*)VT,p,
&zzero,gtightf+rs*p*q*R, p);
}
LTFAT_SAFEFREEALL(gfwork,Sf,S,U,VT);
}
|
