You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
* Tests whether an operation should be applied to the left side.
37
+
*
38
+
* @private
39
+
* @param {string} side - operation side
40
+
* @returns {boolean} boolean indicating if an operation should be applied to the left side
41
+
*/
42
+
functionisLeftSide(side){
43
+
returnside==='left';
44
+
}
45
+
46
+
33
47
// MAIN //
34
48
35
49
/**
36
50
* Applies a real elementary reflector `H = I - tau * v * v^T` to a real M by N matrix `C`.
37
51
*
38
52
* ## Notes
39
53
*
40
-
* - `work` should have `N` indexed elements if side = `'left'` and `M` indexed elements if side = `'right'`.
41
-
* - `V` should have `1 + (M-1) * abs(strideV)` indexed elements if side = `'left'` and `1 + (N-1) * abs(strideV)` indexed elements if side = `'right'`.
42
-
* - `C` is overwritten by `H * C` if side = `'left'` and `C * H` if side = `'right'`.
54
+
* - If `side = 'left'`,
55
+
*
56
+
* - `work` should have `N` indexed elements.
57
+
* - `V` should have `1 + (M-1) * abs(strideV)` indexed elements.
58
+
* - `C` is overwritten by `H * C`.
59
+
*
60
+
* - If `side = 'right'`,
61
+
*
62
+
* - `work` should have `M` indexed elements.
63
+
* - `V` should have `1 + (N-1) * abs(strideV)` indexed elements.
64
+
* - `C` is overwritten by `C * H`.
43
65
*
44
66
* @private
45
-
* @param {string} side - specifies the side of multiplication with `C`. Use `'left'` to form `H * C` and `'right'` to form `C * H`.
67
+
* @param {string} side - specifies the side of multiplication with `C`
46
68
* @param {NonNegativeInteger} M - number of rows in `C`
47
69
* @param {NonNegativeInteger} N - number of columns in `C`
48
70
* @param {Float64Array} V - the vector `v`
@@ -73,51 +95,85 @@ function dlarf1f( side, M, N, V, strideV, offsetV, tau, C, strideC1, strideC2, o
73
95
varlastc;
74
96
vari;
75
97
76
-
lastv=1;
98
+
if(tau===0.0){
99
+
returnC;
100
+
}
77
101
lastc=0;
78
-
if(tau!==0.0){
79
-
if(side==='left'){
80
-
lastv=M;
81
-
}else{
82
-
lastv=N;
83
-
}
84
-
85
-
// i points to the last element in V
86
-
i=offsetV+((lastv-1)*strideV);
102
+
if(isLeftSide(side)){
103
+
lastv=M;
104
+
}else{
105
+
lastv=N;
106
+
}
107
+
// Initialize `i` to point to the last element in `V`:
108
+
i=offsetV+((lastv-1)*strideV);
87
109
88
-
// Move i to the last non-zero element in V
89
-
while(lastv>0&&V[i]===0.0){
90
-
lastv-=1;
91
-
i-=strideV;
92
-
}
93
-
if(side==='left'){
94
-
lastc=iladlc(lastv+1,N,C,strideC1,strideC2,offsetC)+1;// to account for the difference between zero-based and one-based indexing
95
-
}else{
96
-
lastc=iladlr(M,lastv+1,C,strideC1,strideC2,offsetC)+1;// to account for the difference between zero-based and one-based indexing
97
-
}
98
-
// lastc is zero if a matrix is filled with zeros
110
+
// Move `i` to the last non-zero element in `V`, where we assume that V[0] = 1, and it is not stored, so we shouldn't access it...
111
+
while(lastv>0&&V[i]===0.0){
112
+
lastv-=1;
113
+
i-=strideV;
114
+
}
115
+
if(isLeftSide(side)){
116
+
// Scan for the last non-zero column in `C`:
117
+
lastc=iladlc(lastv+1,N,C,strideC1,strideC2,offsetC)+1;// adjust by `+1` to account for the difference between zero-based and one-based indexing
118
+
}else{
119
+
// Scan for the last non-zero row in `C`:
120
+
lastc=iladlr(M,lastv+1,C,strideC1,strideC2,offsetC)+1;// // adjust by `+1` to account for the difference between zero-based and one-based indexing
99
121
}
122
+
// Return `C` unchanged if all elements in `C` are zero...
100
123
if(lastc===0){
101
-
// Returns C unchanged if tau is zero or all elements in C are zero
102
124
returnC;
103
125
}
104
-
if(side==='left'){
126
+
if(isLeftSide(side)){
127
+
// Form: H*C
128
+
129
+
// If `lastv = 1`, this means `V = 1`, so we just need to compute `C = H*C = (1\tau)*C`...
105
130
if(lastv===0){
106
-
dscal(lastc,1.0-tau,C,strideC2,offsetC);// scale the first row
131
+
// C[0,0:lastc] = (1-tau)*C[0,0:lastc]
132
+
dscal(lastc,1.0-tau,C,strideC2,offsetC);// scale the first row
107
133
}else{
108
-
dgemv('transpose',lastv-1,lastc,1.0,C,strideC1,strideC2,offsetC+strideC1,V,strideV,offsetV+strideV,0.0,work,strideWork,offsetWork);// C( 1, 0 ) is accessed here
dscal(lastc,1.0-tau,C,strideC1,offsetC);// scale the first column
115
-
}else{
116
-
dgemv('no-transpose',lastc,lastv-1,1.0,C,strideC1,strideC2,offsetC+strideC2,V,strideV,offsetV+strideV,0.0,work,strideWork,offsetWork);// C( 0, 1 ) is accessed here
117
-
daxpy(lastc,1.0,C,strideC1,offsetC,work,strideWork,offsetWork);// operates on the first column of C
118
-
daxpy(lastc,-tau,work,strideWork,offsetWork,C,strideC1,offsetC);// operates on the first column of C
119
-
dger(lastc,lastv-1,-tau,work,strideWork,offsetWork,V,strideV,offsetV+strideV,C,strideC1,strideC2,offsetC+strideC2);// C( 0, 1 ) is accessed here
150
+
returnC;
151
+
}
152
+
// side === 'right'
153
+
154
+
// Form: C*H
155
+
156
+
// If `N = 1`, then `V = 1`, so we just need to compute `C = CH = C*(1-tau)`...
157
+
if(lastv===0){
158
+
// C[0:lastc,0] = ( 1-tau ) * C[0:lastc,0]
159
+
dscal(lastc,1.0-tau,C,strideC1,offsetC);// scale the first column
@@ -33,11 +36,19 @@ var base = require( './base.js' );
33
36
*
34
37
* ## Notes
35
38
*
36
-
* - `work` should have `N` indexed elements if side = `'left'` and `M` indexed elements if side = `'right'`.
37
-
* - `V` should have `1 + (M-1) * abs(strideV)` indexed elements if side = `'left'` and `1 + (N-1) * abs(strideV)` indexed elements if side = `'right'`.
38
-
* - `C` is overwritten by `H * C` if side = `'left'` and `C * H` if side = `'right'`.
39
+
* - If `side = 'left'`,
39
40
*
40
-
* @param {string} side - specifies the side of multiplication with `C`. Use `'left'` to form `H * C` and `'right'` to form `C * H`.
41
+
* - `work` should have `N` indexed elements.
42
+
* - `V` should have `1 + (M-1) * abs(strideV)` indexed elements.
43
+
* - `C` is overwritten by `H * C`.
44
+
*
45
+
* - If `side = 'right'`,
46
+
*
47
+
* - `work` should have `M` indexed elements.
48
+
* - `V` should have `1 + (N-1) * abs(strideV)` indexed elements.
49
+
* - `C` is overwritten by `C * H`.
50
+
*
51
+
* @param {string} side - specifies the side of multiplication with `C`
41
52
* @param {NonNegativeInteger} M - number of rows in `C`
42
53
* @param {NonNegativeInteger} N - number of columns in `C`
43
54
* @param {Float64Array} V - the vector `v`
@@ -65,8 +76,8 @@ var base = require( './base.js' );
0 commit comments