@@ -329,25 +329,33 @@ end
329
329
330
330
const BiTriSym = Union{Bidiagonal,Tridiagonal,SymTridiagonal}
331
331
const BiTri = Union{Bidiagonal,Tridiagonal}
332
- A_mul_B! (C:: AbstractMatrix , A:: SymTridiagonal , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
333
- A_mul_B! (C:: AbstractMatrix , A:: BiTri , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
334
- A_mul_B! (C:: AbstractMatrix , A:: BiTriSym , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
335
- A_mul_B! (C:: AbstractMatrix , A:: AbstractTriangular , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
336
- A_mul_B! (C:: AbstractMatrix , A:: AbstractMatrix , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
337
- A_mul_B! (C:: AbstractMatrix , A:: Diagonal , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
338
- A_mul_B! (C:: AbstractVector , A:: BiTri , B:: AbstractVector ) = A_mul_B_td! (C, A, B)
339
- A_mul_B! (C:: AbstractMatrix , A:: BiTri , B:: AbstractVecOrMat ) = A_mul_B_td! (C, A, B)
340
- A_mul_B! (C:: AbstractVecOrMat , A:: BiTri , B:: AbstractVecOrMat ) = A_mul_B_td! (C, A, B)
341
-
342
- \ (:: Diagonal , :: RowVector ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
343
- \ (:: Bidiagonal , :: RowVector ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
344
- \ (:: Bidiagonal{<:Number} , :: RowVector{<:Number} ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
345
-
346
- At_ldiv_B (:: Bidiagonal , :: RowVector ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
347
- At_ldiv_B (:: Bidiagonal{<:Number} , :: RowVector{<:Number} ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
348
-
349
- Ac_ldiv_B (:: Bidiagonal , :: RowVector ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
350
- Ac_ldiv_B (:: Bidiagonal{<:Number} , :: RowVector{<:Number} ) = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
332
+ mul! (C:: AbstractMatrix , A:: SymTridiagonal , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
333
+ mul! (C:: AbstractMatrix , A:: BiTri , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
334
+ mul! (C:: AbstractMatrix , A:: BiTriSym , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
335
+ mul! (C:: AbstractMatrix , A:: AbstractTriangular , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
336
+ mul! (C:: AbstractMatrix , A:: AbstractMatrix , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
337
+ mul! (C:: AbstractMatrix , A:: Diagonal , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
338
+ mul! (C:: AbstractMatrix , A:: Adjoint{<:Any,<:Diagonal} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
339
+ mul! (C:: AbstractMatrix , A:: Transpose{<:Any,<:Diagonal} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
340
+ mul! (C:: AbstractMatrix , A:: Adjoint{<:Any,<:AbstractTriangular} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
341
+ mul! (C:: AbstractMatrix , A:: Transpose{<:Any,<:AbstractTriangular} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
342
+ mul! (C:: AbstractMatrix , A:: Adjoint{<:Any,<:AbstractVecOrMat} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
343
+ mul! (C:: AbstractMatrix , A:: Transpose{<:Any,<:AbstractVecOrMat} , B:: BiTriSym ) = A_mul_B_td! (C, A, B)
344
+ mul! (C:: AbstractVector , A:: BiTri , B:: AbstractVector ) = A_mul_B_td! (C, A, B)
345
+ mul! (C:: AbstractMatrix , A:: BiTri , B:: AbstractVecOrMat ) = A_mul_B_td! (C, A, B)
346
+ mul! (C:: AbstractVecOrMat , A:: BiTri , B:: AbstractVecOrMat ) = A_mul_B_td! (C, A, B)
347
+ mul! (C:: AbstractMatrix , A:: BiTri , B:: Transpose{<:Any,<:AbstractVecOrMat} ) = A_mul_B_td! (C, A, B) # around bidiag line 330
348
+ mul! (C:: AbstractMatrix , A:: BiTri , B:: Adjoint{<:Any,<:AbstractVecOrMat} ) = A_mul_B_td! (C, A, B)
349
+ mul! (C:: AbstractVector , A:: BiTri , B:: Transpose{<:Any,<:AbstractVecOrMat} ) = throw (MethodError (mul!, (C, A, B)))
350
+
351
+ \ (:: Diagonal , :: RowVector ) = _mat_ldiv_rowvec_error ()
352
+ \ (:: Bidiagonal , :: RowVector ) = _mat_ldiv_rowvec_error ()
353
+ \ (:: Bidiagonal{<:Number} , :: RowVector{<:Number} ) = _mat_ldiv_rowvec_error ()
354
+ \ (:: Adjoint{<:Any,<:Bidiagonal} , :: RowVector ) = _mat_ldiv_rowvec_error ()
355
+ \ (:: Transpose{<:Any,<:Bidiagonal} , :: RowVector ) = _mat_ldiv_rowvec_error ()
356
+ \ (:: Adjoint{<:Number,<:Bidiagonal{<:Number}} , :: RowVector{<:Number} ) = _mat_ldiv_rowvec_error ()
357
+ \ (:: Transpose{<:Number,<:Bidiagonal{<:Number}} , :: RowVector{<:Number} ) = _mat_ldiv_rowvec_error ()
358
+ _mat_ldiv_rowvec_error () = throw (DimensionMismatch (" Cannot left-divide matrix by transposed vector" ))
351
359
352
360
function check_A_mul_B!_sizes (C, A, B)
353
361
nA, mA = size (A)
379
387
function A_mul_B_td! (C:: AbstractMatrix , A:: BiTriSym , B:: BiTriSym )
380
388
check_A_mul_B!_sizes (C, A, B)
381
389
n = size (A,1 )
382
- n <= 3 && return A_mul_B ! (C, Array (A), Array (B))
390
+ n <= 3 && return mul ! (C, Array (A), Array (B))
383
391
fill! (C, zero (eltype (C)))
384
392
Al = _diag (A, - 1 )
385
393
Ad = _diag (A, 0 )
@@ -438,7 +446,7 @@ function A_mul_B_td!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat)
438
446
if size (C,2 ) != nB
439
447
throw (DimensionMismatch (" A has second dimension $nA , B has $(size (B,2 )) , C has $(size (C,2 )) but all must match" ))
440
448
end
441
- nA <= 3 && return A_mul_B ! (C, Array (A), Array (B))
449
+ nA <= 3 && return mul ! (C, Array (A), Array (B))
442
450
l = _diag (A, - 1 )
443
451
d = _diag (A, 0 )
444
452
u = _diag (A, 1 )
459
467
function A_mul_B_td! (C:: AbstractMatrix , A:: AbstractMatrix , B:: BiTriSym )
460
468
check_A_mul_B!_sizes (C, A, B)
461
469
n = size (A,1 )
462
- n <= 3 && return A_mul_B ! (C, Array (A), Array (B))
470
+ n <= 3 && return mul ! (C, Array (A), Array (B))
463
471
m = size (B,2 )
464
472
Bl = _diag (B, - 1 )
465
473
Bd = _diag (B, 0 )
@@ -493,15 +501,17 @@ const SpecialMatrix = Union{Bidiagonal,SymTridiagonal,Tridiagonal}
493
501
* (A:: SpecialMatrix , B:: SpecialMatrix ) = Array (A) * Array (B)
494
502
495
503
# Generic multiplication
496
- for func in (:* , :Ac_mul_B , :A_mul_Bc , :/ , :A_rdiv_Bc )
497
- @eval ($ func)(A:: Bidiagonal{T} , B:: AbstractVector{T} ) where {T} = ($ func)(Array (A), B)
498
- end
504
+ * (A:: Bidiagonal{T} , B:: AbstractVector{T} ) where {T} = * (Array (A), B)
505
+ * (adjA:: Adjoint{<:Any,<:Bidiagonal{T}} , B:: AbstractVector{T} ) where {T} = * (Adjoint (Array (adjA. parent)), B)
506
+ * (A:: Bidiagonal{T} , adjB:: Adjoint{<:Any,<:AbstractVector{T}} ) where {T} = * (Array (A), Adjoint (adjB. parent))
507
+ / (A:: Bidiagonal{T} , B:: AbstractVector{T} ) where {T} = / (Array (A), B)
508
+ / (A:: Bidiagonal{T} , adjB:: Adjoint{<:Any,<:AbstractVector{T}} ) where {T} = / (Array (A), Adjoint (adjB. parent))
499
509
500
510
# Linear solvers
501
- A_ldiv_B ! (A:: Union{Bidiagonal, AbstractTriangular} , b:: AbstractVector ) = naivesub! (A, b)
502
- At_ldiv_B! (A :: Bidiagonal , b:: AbstractVector ) = A_ldiv_B ! (transpose (A ), b)
503
- Ac_ldiv_B! (A :: Bidiagonal , b:: AbstractVector ) = A_ldiv_B ! (adjoint (A ), b)
504
- function A_ldiv_B ! (A:: Union{Bidiagonal,AbstractTriangular} , B:: AbstractMatrix )
511
+ ldiv ! (A:: Union{Bidiagonal, AbstractTriangular} , b:: AbstractVector ) = naivesub! (A, b)
512
+ ldiv! (transA :: Transpose{<:Any,<: Bidiagonal} , b:: AbstractVector ) = ldiv ! (transpose (transA . parent ), b)
513
+ ldiv! (adjA :: Adjoint{<:Any,<: Bidiagonal} , b:: AbstractVector ) = ldiv ! (adjoint (adjA . parent ), b)
514
+ function ldiv ! (A:: Union{Bidiagonal,AbstractTriangular} , B:: AbstractMatrix )
505
515
nA,mA = size (A)
506
516
tmp = similar (B,size (B,1 ))
507
517
n = size (B, 1 )
@@ -510,26 +520,40 @@ function A_ldiv_B!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
510
520
end
511
521
for i = 1 : size (B,2 )
512
522
copy! (tmp, 1 , B, (i - 1 )* n + 1 , n)
513
- A_ldiv_B ! (A, tmp)
523
+ ldiv ! (A, tmp)
514
524
copy! (B, (i - 1 )* n + 1 , tmp, 1 , n) # Modify this when array view are implemented.
515
525
end
516
526
B
517
527
end
518
- for func in (:Ac_ldiv_B! , :At_ldiv_B! )
519
- @eval function ($ func)(A:: Union{Bidiagonal,AbstractTriangular} , B:: AbstractMatrix )
520
- nA,mA = size (A)
521
- tmp = similar (B,size (B,1 ))
522
- n = size (B, 1 )
523
- if mA != n
524
- throw (DimensionMismatch (" size of A' is ($mA ,$nA ), corresponding dimension of B is $n " ))
525
- end
526
- for i = 1 : size (B,2 )
527
- copy! (tmp, 1 , B, (i - 1 )* n + 1 , n)
528
- ($ func)(A, tmp)
529
- copy! (B, (i - 1 )* n + 1 , tmp, 1 , n) # Modify this when array view are implemented.
530
- end
531
- B
528
+ function ldiv! (adjA:: Adjoint{<:Any,<:Union{Bidiagonal,AbstractTriangular}} , B:: AbstractMatrix )
529
+ A = adjA. parent
530
+ nA,mA = size (A)
531
+ tmp = similar (B,size (B,1 ))
532
+ n = size (B, 1 )
533
+ if mA != n
534
+ throw (DimensionMismatch (" size of A' is ($mA ,$nA ), corresponding dimension of B is $n " ))
532
535
end
536
+ for i = 1 : size (B,2 )
537
+ copy! (tmp, 1 , B, (i - 1 )* n + 1 , n)
538
+ ldiv! (Adjoint (A), tmp)
539
+ copy! (B, (i - 1 )* n + 1 , tmp, 1 , n) # Modify this when array view are implemented.
540
+ end
541
+ B
542
+ end
543
+ function ldiv! (transA:: Transpose{<:Any,<:Union{Bidiagonal,AbstractTriangular}} , B:: AbstractMatrix )
544
+ A = transA. parent
545
+ nA,mA = size (A)
546
+ tmp = similar (B,size (B,1 ))
547
+ n = size (B, 1 )
548
+ if mA != n
549
+ throw (DimensionMismatch (" size of A' is ($mA ,$nA ), corresponding dimension of B is $n " ))
550
+ end
551
+ for i = 1 : size (B,2 )
552
+ copy! (tmp, 1 , B, (i - 1 )* n + 1 , n)
553
+ ldiv! (Transpose (A), tmp)
554
+ copy! (B, (i - 1 )* n + 1 , tmp, 1 , n) # Modify this when array view are implemented.
555
+ end
556
+ B
533
557
end
534
558
# Generic solver using naive substitution
535
559
function naivesub! (A:: Bidiagonal{T} , b:: AbstractVector , x:: AbstractVector = b) where T
@@ -554,15 +578,26 @@ function naivesub!(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b) w
554
578
end
555
579
556
580
# ## Generic promotion methods and fallbacks
557
- for (f,g) in ((:\ , :A_ldiv_B! ), (:At_ldiv_B , :At_ldiv_B! ), (:Ac_ldiv_B , :Ac_ldiv_B! ))
558
- @eval begin
559
- function ($ f)(A:: Bidiagonal{TA} , B:: AbstractVecOrMat{TB} ) where {TA<: Number ,TB<: Number }
560
- TAB = typeof ((zero (TA)* zero (TB) + zero (TA)* zero (TB))/ one (TA))
561
- ($ g)(convert (AbstractArray{TAB}, A), copy_oftype (B, TAB))
562
- end
563
- ($ f)(A:: Bidiagonal , B:: AbstractVecOrMat ) = ($ g)(A, copy (B))
564
- end
581
+ function \ (A:: Bidiagonal{<:Number} , B:: AbstractVecOrMat{<:Number} )
582
+ TA, TB = eltype (A), eltype (B)
583
+ TAB = typeof ((zero (TA)* zero (TB) + zero (TA)* zero (TB))/ one (TA))
584
+ ldiv! (convert (AbstractArray{TAB}, A), copy_oftype (B, TAB))
585
+ end
586
+ \ (A:: Bidiagonal , B:: AbstractVecOrMat ) = ldiv! (A, copy (B))
587
+ function \ (transA:: Transpose{<:Number,<:Bidiagonal{<:Number}} , B:: AbstractVecOrMat{<:Number} )
588
+ A = transA. parent
589
+ TA, TB = eltype (A), eltype (B)
590
+ TAB = typeof ((zero (TA)* zero (TB) + zero (TA)* zero (TB))/ one (TA))
591
+ ldiv! (Transpose (convert (AbstractArray{TAB}, A)), copy_oftype (B, TAB))
592
+ end
593
+ \ (transA:: Transpose{<:Any,<:Bidiagonal} , B:: AbstractVecOrMat ) = ldiv! (Transpose (transA. parent), copy (B))
594
+ function \ (adjA:: Adjoint{<:Number,<:Bidiagonal{<:Number}} , B:: AbstractVecOrMat{<:Number} )
595
+ A = adjA. parent
596
+ TA, TB = eltype (A), eltype (B)
597
+ TAB = typeof ((zero (TA)* zero (TB) + zero (TA)* zero (TB))/ one (TA))
598
+ ldiv! (Adjoint (convert (AbstractArray{TAB}, A)), copy_oftype (B, TAB))
565
599
end
600
+ \ (adjA:: Adjoint{<:Any,<:Bidiagonal} , B:: AbstractVecOrMat ) = ldiv! (Adjoint (adjA. parent), copy (B))
566
601
567
602
factorize (A:: Bidiagonal ) = A
568
603
0 commit comments