|
63 | 63 | """
|
64 | 64 | function ProductFun(cfs::AbstractMatrix{T},sp::AbstractProductSpace{Tuple{S,V},DD};
|
65 | 65 | tol::Real=100eps(T),chopping::Bool=false) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,DD}
|
| 66 | + |
| 67 | + kend = size(cfs, 2) |
66 | 68 | if chopping
|
67 |
| - ncfs, kend = norm(cfs,Inf), size(cfs,2) |
68 |
| - if kend > 1 |
69 |
| - while kend > 0 |
70 |
| - if all(iszero, @view(cfs[:, kend])) |
71 |
| - kend -= 1 |
72 |
| - continue |
73 |
| - end |
74 |
| - if !isempty(chop(@view(cfs[:,kend]), ncfs*tol)) |
75 |
| - break |
76 |
| - end |
77 |
| - kend-=1 |
78 |
| - end |
79 |
| - end |
80 |
| - ret=VFun{S,T}[Fun(columnspace(sp,k),chop(@view(cfs[:,k]), ncfs*tol)) for k=1:max(kend,1)] |
81 |
| - ProductFun{S,V,typeof(sp),T}(ret,sp) |
| 69 | + ncfs = norm(cfs,Inf) |
| 70 | + kend -= ntrailingzerocols(cfs, ncfs*tol) |
| 71 | + end |
| 72 | + |
| 73 | + ret = if kend == 0 |
| 74 | + VFun{S,T}[Fun(columnspace(sp,1), T[]) for k=1:1] |
82 | 75 | else
|
83 |
| - ret=VFun{S,T}[Fun(columnspace(sp,k),cfs[:,k]) for k=1:size(cfs,2)] |
84 |
| - ProductFun{S,V,typeof(sp),T}(ret,sp) |
| 76 | + VFun{S,T}[Fun(columnspace(sp,k), |
| 77 | + if chopping |
| 78 | + chop(@view(cfs[:,k]), ncfs*tol) |
| 79 | + else |
| 80 | + cfs[:,k] |
| 81 | + end |
| 82 | + ) |
| 83 | + for k=1:kend |
| 84 | + ] |
85 | 85 | end
|
| 86 | + |
| 87 | + ProductFun{S,V,typeof(sp),T}(ret,sp) |
86 | 88 | end
|
87 | 89 |
|
88 | 90 | ## Construction in a ProductSpace via a Vector of Funs
|
@@ -147,7 +149,32 @@ ProductFun(f::Function) = ProductFun(dynamic(f),ChebyshevInterval(),ChebyshevInt
|
147 | 149 | ## Conversion from other 2D Funs
|
148 | 150 |
|
149 | 151 | ProductFun(f::LowRankFun; kw...) = ProductFun(coefficients(f),space(f,1),space(f,2); kw...)
|
150 |
| -ProductFun(f::Fun{S}; kw...) where {S<:AbstractProductSpace} = ProductFun(coefficientmatrix(f), space(f); kw...) |
| 152 | +function nzerofirst(itr, tol=0.0) |
| 153 | + n = 0 |
| 154 | + for v in itr |
| 155 | + if all(iszero, v) || (tol == 0 ? false : isempty(chop(v, tol))) |
| 156 | + n += 1 |
| 157 | + else |
| 158 | + break |
| 159 | + end |
| 160 | + end |
| 161 | + n |
| 162 | +end |
| 163 | +function ntrailingzerocols(A, tol=0.0) |
| 164 | + itr = (view(A, :, i) for i in reverse(axes(A,2))) |
| 165 | + nzerofirst(itr, tol) |
| 166 | +end |
| 167 | +function ntrailingzerorows(A, tol=0.0) |
| 168 | + itr = (view(A, i, :) for i in reverse(axes(A,1))) |
| 169 | + nzerofirst(itr, tol) |
| 170 | +end |
| 171 | +function ProductFun(f::Fun{<:AbstractProductSpace}; kw...) |
| 172 | + M = coefficientmatrix(f) |
| 173 | + nc = ntrailingzerocols(M) |
| 174 | + nr = ntrailingzerorows(M) |
| 175 | + A = @view M[1:max(1,end-nr), 1:max(1,end-nc)] |
| 176 | + ProductFun(A, space(f); kw...) |
| 177 | +end |
151 | 178 |
|
152 | 179 | ## Conversion to other ProductSpaces with the same coefficients
|
153 | 180 |
|
@@ -181,7 +208,7 @@ ProductFun(f::Fun,sp::BivariateSpace) = ProductFun([Fun(f,columnspace(sp,1))],sp
|
181 | 208 |
|
182 | 209 |
|
183 | 210 | function funlist2coefficients(f::Vector{VFun{S,T}}) where {S,T}
|
184 |
| - A=zeros(T,mapreduce(ncoefficients,max,f),length(f)) |
| 211 | + A=zeros(T,mapreduce(ncoefficients,max,f,init=0),length(f)) |
185 | 212 | for k=1:length(f)
|
186 | 213 | A[1:ncoefficients(f[k]),k]=f[k].coefficients
|
187 | 214 | end
|
|
0 commit comments