-
-
Notifications
You must be signed in to change notification settings - Fork 46
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Refactor linear algebra module and remove unused code #195
base: main
Are you sure you want to change the base?
Changes from 1 commit
6f4fe35
d667c95
10eea19
5bd805c
f33a31b
3a0af32
11ed075
c3d1e2e
ea376ac
c7356d7
0832f6d
ae0628b
2a2e7d8
e82bda8
88776da
b84f2ba
ab34e62
14d3f67
02da167
b2f481c
e43bfa5
bfa6907
961475e
f6cb782
569a96b
d8a4fc2
c66401f
d399e1a
1d5e441
9f46519
d376d36
99a3a2b
5f9a1fe
e83de29
b06e436
36ae80e
63adcf6
8afd14a
3af89bf
d9fd254
79d8a8d
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
…ormat
- Loading branch information
There are no files selected for viewing
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,70 @@ | ||
module lapack | ||
|
||
import vsl.lapack.lapack64 | ||
|
||
// Direct specifies the direction of the multiplication for the Householder matrix. | ||
pub type Direct = lapack64.Direct | ||
|
||
// Sort is the sorting order. | ||
pub type Sort = lapack64.Sort | ||
|
||
// StoreV indicates the storage direction of elementary reflectors. | ||
pub type StoreV = lapack64.StoreV | ||
|
||
// MatrixNorm represents the kind of matrix norm to compute. | ||
pub type MatrixNorm = lapack64.MatrixNorm | ||
|
||
// MatrixType represents the kind of matrix represented in the data. | ||
pub type MatrixType = lapack64.MatrixType | ||
|
||
// Pivot specifies the pivot type for plane rotations. | ||
pub type Pivot = lapack64.Pivot | ||
|
||
// ApplyOrtho specifies which orthogonal matrix is applied in Dormbr. | ||
pub type ApplyOrtho = lapack64.ApplyOrtho | ||
|
||
// GenOrtho specifies which orthogonal matrix is generated in Dorgbr. | ||
pub type GenOrtho = lapack64.GenOrtho | ||
|
||
// SVDJob specifies the singular vector computation type for SVD. | ||
pub type SVDJob = lapack64.SVDJob | ||
|
||
// GSVDJob specifies the singular vector computation type for Generalized SVD. | ||
pub type GSVDJob = lapack64.GSVDJob | ||
|
||
// EVComp specifies how eigenvectors are computed in Dsteqr. | ||
pub type EVComp = lapack64.EVComp | ||
|
||
// EVJob specifies whether eigenvectors are computed in Dsyev. | ||
pub type EVJob = lapack64.EVJob | ||
|
||
// LeftEVJob specifies whether left eigenvectors are computed in Dgeev. | ||
pub type LeftEVJob = lapack64.LeftEVJob | ||
|
||
// RightEVJob specifies whether right eigenvectors are computed in Dgeev. | ||
pub type RightEVJob = lapack64.RightEVJob | ||
|
||
// BalanceJob specifies matrix balancing operation. | ||
pub type BalanceJob = lapack64.BalanceJob | ||
|
||
// SchurJob specifies whether the Schur form is computed in Dhseqr. | ||
pub type SchurJob = lapack64.SchurJob | ||
|
||
// SchurComp specifies whether and how the Schur vectors are computed in Dhseqr. | ||
pub type SchurComp = lapack64.SchurComp | ||
|
||
// UpdateSchurComp specifies whether the matrix of Schur vectors is updated in Dtrexc. | ||
pub type UpdateSchurComp = lapack64.UpdateSchurComp | ||
|
||
// EVSide specifies what eigenvectors are computed in Dtrevc3. | ||
pub type EVSide = lapack64.EVSide | ||
|
||
// EVHowMany specifies which eigenvectors are computed in Dtrevc3 and how. | ||
pub type EVHowMany = lapack64.EVHowMany | ||
|
||
// MaximizeNormXJob specifies the heuristic method for computing a contribution to | ||
// the reciprocal Dif-estimate in Dlatdf. | ||
pub type MaximizeNormXJob = lapack64.MaximizeNormXJob | ||
|
||
// OrthoComp specifies whether and how the orthogonal matrix is computed in Dgghrd. | ||
pub type OrthoComp = lapack64.OrthoComp |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,199 @@ | ||
module lapack64 | ||
|
||
// Direct specifies the direction of the multiplication for the Householder matrix. | ||
pub enum Direct as u8 { | ||
// Reflectors are right-multiplied, H_0 * H_1 * ... * H_{k-1}. | ||
forward = u8(`F`) | ||
// Reflectors are left-multiplied, H_{k-1} * ... * H_1 * H_0. | ||
backward = u8(`B`) | ||
} | ||
|
||
// Sort is the sorting order. | ||
pub enum Sort as u8 { | ||
sort_increasing = u8(`I`) | ||
sort_decreasing = u8(`D`) | ||
} | ||
|
||
// StoreV indicates the storage direction of elementary reflectors. | ||
pub enum StoreV as u8 { | ||
// Reflector stored in a column of the matrix. | ||
column_wise = u8(`C`) | ||
// Reflector stored in a row of the matrix. | ||
row_wise = u8(`R`) | ||
} | ||
|
||
// MatrixNorm represents the kind of matrix norm to compute. | ||
pub enum MatrixNorm as u8 { | ||
// max(abs(A(i,j))) | ||
max_abs = u8(`M`) | ||
// Maximum absolute column sum (one norm) | ||
max_column_sum = u8(`O`) | ||
// Maximum absolute row sum (infinity norm) | ||
max_row_sum = u8(`I`) | ||
// Frobenius norm (sqrt of sum of squares) | ||
frobenius = u8(`F`) | ||
} | ||
|
||
// MatrixType represents the kind of matrix represented in the data. | ||
pub enum MatrixType as u8 { | ||
// A general dense matrix. | ||
general = u8(`G`) | ||
// An upper triangular matrix. | ||
upper_tri = u8(`U`) | ||
// A lower triangular matrix. | ||
lower_tri = u8(`L`) | ||
} | ||
|
||
// Pivot specifies the pivot type for plane rotations. | ||
pub enum Pivot as u8 { | ||
variable = u8(`V`) | ||
top = u8(`T`) | ||
bottom = u8(`B`) | ||
} | ||
|
||
// ApplyOrtho specifies which orthogonal matrix is applied in Dormbr. | ||
pub enum ApplyOrtho as u8 { | ||
// Apply P or Pᵀ. | ||
apply_p = u8(`P`) | ||
// Apply Q or Qᵀ. | ||
apply_q = u8(`Q`) | ||
} | ||
|
||
// GenOrtho specifies which orthogonal matrix is generated in Dorgbr. | ||
pub enum GenOrtho as u8 { | ||
// Generate Pᵀ. | ||
generate_pt = u8(`P`) | ||
// Generate Q. | ||
generate_q = u8(`Q`) | ||
} | ||
|
||
// SVDJob specifies the singular vector computation type for SVD. | ||
pub enum SVDJob as u8 { | ||
// Compute all columns of the orthogonal matrix U or V. | ||
svd_all = u8(`A`) | ||
// Compute the singular vectors and store them in the orthogonal matrix U or V. | ||
svd_store = u8(`S`) | ||
// Compute the singular vectors and overwrite them on the input matrix A. | ||
svd_overwrite = u8(`O`) | ||
// Do not compute singular vectors. | ||
svd_none = u8(`N`) | ||
} | ||
|
||
// GSVDJob specifies the singular vector computation type for Generalized SVD. | ||
pub enum GSVDJob as u8 { | ||
// Compute orthogonal matrix U. | ||
gsvd_u = u8(`U`) | ||
// Compute orthogonal matrix V. | ||
gsvd_v = u8(`V`) | ||
// Compute orthogonal matrix Q. | ||
gsvd_q = u8(`Q`) | ||
// Use unit-initialized matrix. | ||
gsvd_unit = u8(`I`) | ||
// Do not compute orthogonal matrix. | ||
gsvd_none = u8(`N`) | ||
} | ||
|
||
// EVComp specifies how eigenvectors are computed in Dsteqr. | ||
pub enum EVComp as u8 { | ||
// Compute eigenvectors of the original symmetric matrix. | ||
ev_orig = u8(`V`) | ||
// Compute eigenvectors of the tridiagonal matrix. | ||
ev_tridiag = u8(`I`) | ||
// Do not compute eigenvectors. | ||
ev_comp_none = u8(`N`) | ||
} | ||
|
||
// EVJob specifies whether eigenvectors are computed in Dsyev. | ||
pub enum EVJob as u8 { | ||
// Compute eigenvectors. | ||
ev_compute = u8(`V`) | ||
// Do not compute eigenvectors. | ||
ev_none = u8(`N`) | ||
} | ||
|
||
// LeftEVJob specifies whether left eigenvectors are computed in Dgeev. | ||
pub enum LeftEVJob as u8 { | ||
// Compute left eigenvectors. | ||
left_ev_compute = u8(`V`) | ||
// Do not compute left eigenvectors. | ||
left_ev_none = u8(`N`) | ||
} | ||
|
||
// RightEVJob specifies whether right eigenvectors are computed in Dgeev. | ||
pub enum RightEVJob as u8 { | ||
// Compute right eigenvectors. | ||
right_ev_compute = u8(`V`) | ||
// Do not compute right eigenvectors. | ||
right_ev_none = u8(`N`) | ||
} | ||
|
||
// BalanceJob specifies matrix balancing operation. | ||
pub enum BalanceJob as u8 { | ||
permute = u8(`P`) | ||
scale = u8(`S`) | ||
permute_scale = u8(`B`) | ||
balance_none = u8(`N`) | ||
} | ||
|
||
// SchurJob specifies whether the Schur form is computed in Dhseqr. | ||
pub enum SchurJob as u8 { | ||
eigenvalues_only = u8(`E`) | ||
eigenvalues_and_schur = u8(`S`) | ||
} | ||
|
||
// SchurComp specifies whether and how the Schur vectors are computed in Dhseqr. | ||
pub enum SchurComp as u8 { | ||
// Compute Schur vectors of the original matrix. | ||
schur_orig = u8(`V`) | ||
// Compute Schur vectors of the upper Hessenberg matrix. | ||
schur_hess = u8(`I`) | ||
// Do not compute Schur vectors. | ||
schur_none = u8(`N`) | ||
} | ||
|
||
// UpdateSchurComp specifies whether the matrix of Schur vectors is updated in Dtrexc. | ||
pub enum UpdateSchurComp as u8 { | ||
// Update the matrix of Schur vectors. | ||
update_schur = u8(`V`) | ||
// Do not update the matrix of Schur vectors. | ||
update_schur_none = u8(`N`) | ||
} | ||
|
||
// EVSide specifies what eigenvectors are computed in Dtrevc3. | ||
pub enum EVSide as u8 { | ||
// Compute only right eigenvectors. | ||
ev_right = u8(`R`) | ||
// Compute only left eigenvectors. | ||
ev_left = u8(`L`) | ||
// Compute both right and left eigenvectors. | ||
ev_both = u8(`B`) | ||
} | ||
|
||
// EVHowMany specifies which eigenvectors are computed in Dtrevc3 and how. | ||
pub enum EVHowMany as u8 { | ||
// Compute all right and/or left eigenvectors. | ||
ev_all = u8(`A`) | ||
// Compute all right and/or left eigenvectors multiplied by an input matrix. | ||
ev_all_mul_q = u8(`B`) | ||
// Compute selected right and/or left eigenvectors. | ||
ev_selected = u8(`S`) | ||
} | ||
|
||
// MaximizeNormXJob specifies the heuristic method for computing a contribution to | ||
// the reciprocal Dif-estimate in Dlatdf. | ||
pub enum MaximizeNormXJob as u8 { | ||
// Solve Z*x=h-f where h is a vector of ±1. | ||
local_look_ahead = 0 | ||
// Compute an approximate null-vector e of Z, normalize e and solve Z*x=±e-f. | ||
normalized_null_vector = 2 | ||
} | ||
|
||
// OrthoComp specifies whether and how the orthogonal matrix is computed in Dgghrd. | ||
pub enum OrthoComp as u8 { | ||
// Do not compute the orthogonal matrix. | ||
ortho_none = u8(`N`) | ||
// The orthogonal matrix is formed explicitly and returned in the argument. | ||
ortho_explicit = u8(`I`) | ||
// The orthogonal matrix is post-multiplied into the matrix stored in the argument on entry. | ||
ortho_postmul = u8(`V`) | ||
} | ||
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,33 @@ | ||
module lapack64 | ||
|
||
import math | ||
import vsl.blas | ||
|
||
// dgebal balances a general real matrix A. | ||
pub fn dgebal(job BalanceJob, n int, mut a []f64, lda int, scale []f64) int { | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
mut info := 0 | ||
if job != .balance_none && job != .permute && job != .scale && job != .permute_scale { | ||
info = -1 | ||
} else if n < 0 { | ||
info = -2 | ||
} else if lda < math.max(1, n) { | ||
info = -4 | ||
} | ||
|
||
if info != 0 { | ||
return info | ||
} | ||
|
||
// Quick return if possible | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
// Placeholder for the actual LAPACK function calls | ||
// Example: info = dgebal(job, n, a, lda, scale) | ||
return info | ||
Comment on lines
+7
to
+32
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Ensure the implementation of the LAPACK function calls is completed as indicated by the placeholder comment. The function currently does not perform matrix balancing but returns an uninitialized |
||
} |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,39 @@ | ||
module lapack64 | ||
|
||
import math | ||
import vsl.blas | ||
|
||
// dgeev computes the eigenvalues and, optionally, the left and/or right eigenvectors for a real nonsymmetric matrix A. | ||
pub fn dgeev(jobvl LeftEVJob, jobvr LeftEVJob, n int, mut a []f64, lda int, wr []f64, wi []f64, mut vl []f64, ldvl int, mut vr []f64, ldvr int) int { | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
mut info := 0 | ||
if jobvl != .left_ev_none && jobvl != .left_ev_compute { | ||
info = -1 | ||
} else if jobvr != .left_ev_none && jobvr != .left_ev_compute { | ||
info = -2 | ||
} else if n < 0 { | ||
info = -3 | ||
} else if lda < math.max(1, n) { | ||
info = -5 | ||
} else if ldvl < 1 || (jobvl == .left_ev_compute && ldvl < n) { | ||
info = -8 | ||
} else if ldvr < 1 || (jobvr == .left_ev_compute && ldvr < n) { | ||
info = -10 | ||
} | ||
|
||
if info != 0 { | ||
return info | ||
} | ||
|
||
// Quick return if possible | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
// Placeholder for the actual LAPACK function calls | ||
// Example: info = dgehrd(n, ilo, ihi, a, lda, tau, work, lwork) | ||
return info | ||
} |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,35 @@ | ||
module lapack64 | ||
|
||
import math | ||
import vsl.blas | ||
|
||
// dgehrd reduces a general real matrix A to upper Hessenberg form H by an orthogonal similarity transformation. | ||
pub fn dgehrd(n int, ilo int, ihi int, mut a []f64, lda int, tau []f64) int { | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
mut info := 0 | ||
if n < 0 { | ||
info = -1 | ||
} else if ilo < 1 || ilo > math.max(1, n) { | ||
info = -2 | ||
} else if ihi < math.min(ilo, n) || ihi > n { | ||
info = -3 | ||
} else if lda < math.max(1, n) { | ||
info = -5 | ||
} | ||
|
||
if info != 0 { | ||
return info | ||
} | ||
|
||
// Quick return if possible | ||
if n == 0 { | ||
return 0 | ||
} | ||
|
||
// Placeholder for the actual LAPACK function calls | ||
// Example: info = dgehrd(n, ilo, ihi, a, lda, tau, work, lwork) | ||
return info | ||
Comment on lines
+7
to
+34
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Ensure the implementation of the LAPACK function calls is completed as indicated by the placeholder comment. The function currently does not perform the reduction to Hessenberg form but returns an uninitialized |
||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The enums
MatrixType
,Pivot
,ApplyOrtho
, etc., are comprehensive and cover a wide range of functionalities. Consider adding unit tests to verify the correctness of each enum value's behavior in functions where they are used.Would you like me to help by creating some unit tests for these enums?