wannyk · October 30, 2015 08:04 · miiabi · Sep 15, 2023
diff --git a/SVD.swift b/SVD.swift
 //
 //  SVD.swift
 //
 //  Created by Dongwan Kim on 2015. 10. 30..
 //

 /* MATLAB svd(X)

 For the matrix

    X =
        1    2
        3    4
        5    6
        7    8

 the statement

    [U,S,V] = svd(X)

 produces

    U =
        -0.1525   -0.8226   -0.3945   -0.3800
        -0.3499   -0.4214    0.2428    0.8007
        -0.5474   -0.0201    0.6979   -0.4614
        -0.7448    0.3812   -0.5462    0.0407

    S =
        14.2691         0
              0    0.6268
              0         0
              0         0

    V =
        -0.6414     0.7672
        -0.7672    -0.6414
 */

 import Foundation
 import Accelerate

 func testSVD() {
    let x = [1.0, 3.0, 5.0, 7.0, 2.0, 4.0, 6.0, 8.0]
    let m = 4
    let n = 2
    var u : [Double]
    var s : [Double]
    var v : [Double]
    (u, s, v) = svd(x, m: m, n: n)
    print(u)
    print(s)
    print(v)
 }

 func svd(x:[Double], m:Int, n:Int) -> (u:[Double], s:[Double], v:[Double]) {
    var JOBZ = Int8(UnicodeScalar("A").value)
    var M = __CLPK_integer(m)
    var N = __CLPK_integer(n)
    var A = x
    var LDA = __CLPK_integer(m)
    var S = [__CLPK_doublereal](count: min(m,n), repeatedValue: 0.0)
    var U = [__CLPK_doublereal](count: m*m, repeatedValue: 0.0)
    var LDU = __CLPK_integer(m)
    var VT = [__CLPK_doublereal](count: n*n, repeatedValue: 0.0)
    var LDVT = __CLPK_integer(n)
    let lwork = min(m,n)*(6+4*min(m,n))+max(m,n)
    var WORK = [__CLPK_doublereal](count: lwork, repeatedValue: 0.0)
    var LWORK = __CLPK_integer(lwork)
    var IWORK = [__CLPK_integer](count: 8*min(m,n), repeatedValue: 0)
    var INFO = __CLPK_integer(0)
    dgesdd_(&JOBZ, &M, &N, &A, &LDA, &S, &U, &LDU, &VT, &LDVT, &WORK, &LWORK, &IWORK, &INFO)
    var s = [Double](count: m*n, repeatedValue: 0.0)
    for ni in 0...n-1 {
        s[ni*m+ni] = S[ni]
    }
    var v = [Double](count: n*n, repeatedValue: 0.0)
    vDSP_mtransD(VT, 1, &v, 1, vDSP_Length(n), vDSP_Length(n))
    return (U, s, v)
 }
	//
	// SVD.swift
	//
	// Created by Dongwan Kim on 2015. 10. 30..
	//

	/* MATLAB svd(X)

	For the matrix

	X =
	1 2
	3 4
	5 6
	7 8

	the statement

	[U,S,V] = svd(X)

	produces

	U =
	-0.1525 -0.8226 -0.3945 -0.3800
	-0.3499 -0.4214 0.2428 0.8007
	-0.5474 -0.0201 0.6979 -0.4614
	-0.7448 0.3812 -0.5462 0.0407

	S =
	14.2691 0
	0 0.6268
	0 0
	0 0

	V =
	-0.6414 0.7672
	-0.7672 -0.6414
	*/

	import Foundation
	import Accelerate

	func testSVD() {
	let x = [1.0, 3.0, 5.0, 7.0, 2.0, 4.0, 6.0, 8.0]
	let m = 4
	let n = 2
	var u : [Double]
	var s : [Double]
	var v : [Double]
	(u, s, v) = svd(x, m: m, n: n)
	print(u)
	print(s)
	print(v)
	}

	func svd(x:[Double], m:Int, n:Int) -> (u:[Double], s:[Double], v:[Double]) {
	var JOBZ = Int8(UnicodeScalar("A").value)
	var M = __CLPK_integer(m)
	var N = __CLPK_integer(n)
	var A = x
	var LDA = __CLPK_integer(m)
	var S = [__CLPK_doublereal](count: min(m,n), repeatedValue: 0.0)
	var U = [__CLPK_doublereal](count: m*m, repeatedValue: 0.0)
	var LDU = __CLPK_integer(m)
	var VT = [__CLPK_doublereal](count: n*n, repeatedValue: 0.0)
	var LDVT = __CLPK_integer(n)
	let lwork = min(m,n)(6+4min(m,n))+max(m,n)
	var WORK = [__CLPK_doublereal](count: lwork, repeatedValue: 0.0)
	var LWORK = __CLPK_integer(lwork)
	var IWORK = [__CLPK_integer](count: 8*min(m,n), repeatedValue: 0)
	var INFO = __CLPK_integer(0)
	dgesdd_(&JOBZ, &M, &N, &A, &LDA, &S, &U, &LDU, &VT, &LDVT, &WORK, &LWORK, &IWORK, &INFO)
	var s = [Double](count: m*n, repeatedValue: 0.0)
	for ni in 0...n-1 {
	s[ni*m+ni] = S[ni]
	}
	var v = [Double](count: n*n, repeatedValue: 0.0)
	vDSP_mtransD(VT, 1, &v, 1, vDSP_Length(n), vDSP_Length(n))
	return (U, s, v)
	}