/
linAlg.php
116 lines (109 loc) · 3.02 KB
/
linAlg.php
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
<?php
declare(strict_types=1);
namespace Np\linAlgb;
use Np\core\{blas,lapack};
use Np\matrix, Np\vector;
/**
* Linear Algebra
*
*
*
* @package Np
* @category Scientific Computing
* @author ghost (Shubham Chaudhary)
* @email ghost.jat@gmail.com
* @copyright (c) 2020-2021, Shubham Chaudhary
*/
trait linAlg {
/**
*
* get dot product of m.m | m.v | v.v
*
* @param \Np\matrix|\Np\vector $d
* @return matrix|vector|float
*/
public function dot(matrix|vector $d): matrix|vector|float {
if ($this instanceof matrix) {
if ($d instanceof matrix) {
return $this->dotMatrix($d);
}
return $this->dotVector($d);
}
return blas::dot($this, $d);
}
/**
* get matrix & matrix dot product
* @param \Np\matrix $matrix
* @return \Np\matrix
*/
protected function dotMatrix(matrix $matrix): matrix {
if ($this->checkDimensions($this,$matrix)) {
$ar = self::factory($this->row, $matrix->col);
blas::gemm($this, $matrix, $ar);
return $ar;
}
}
/**
* get dot product of matrix & a vector
* @param \Np\vector $vector
* @return \Np\vector
*/
protected function dotVector(vector $vector): vector {
if ($this->checkDimensions($vector, $this)) {
$mvr = vector::factory($this->col);
blas::gemv($this, $vector, $mvr);
return $mvr;
}
}
/**
*
* Compute the multiplicative inverse of the matrix.
* @return matrix|null
*/
public function inverse(): matrix|null {
if ($this->isSquare()) {
$imat = $this->copy();
$ipiv = vector::factory($this->row, vector::INT);
$lp = lapack::getrf($imat, $ipiv);
if ($lp != 0) {
return null;
}
$lp = lapack::getri($imat, $ipiv);
if ($lp != 0) {
return null;
}
unset($ipiv);
unset($lp);
return $imat;
}
self::_err('Error::invalid Size of matrix!');
}
/**
* FIXEME:-Bug noticed on 10/06/21
* Compute the (Moore-Penrose) pseudo inverse of the general matrix.
* @return matrix|null
*/
public function pseudoInverse(): matrix|null {
$k = min($this->row, $this->col);
$s = vector::factory($k);
$u = self::factory($this->row, $this->row);
$vt = self::factory($this->col, $this->col);
$imat = $this->copy();
$lp = lapack::gesdd($imat, $s, $u, $vt);
if ($lp != 0) {
return null;
}
for ($i = 0; $i < $k; ++$i) {
blas::scale(1.0 / $s->data[$i], $vt->rowAsVector($i));
}
unset($imat);
unset($k);
unset($lp);
unset($s);
$mr = self::factory($this->col, $this->row, $this->dtype);
blas::gemm($vt, $u, $mr);
unset($u);
unset($vt);
return $mr;
}
}