Skip to content

Commit

Permalink
1. added stripping of the real or imaginary parts of the matrix/array,
Browse files Browse the repository at this point in the history
2. finished error checks of sgemm
  • Loading branch information
jacobbogers committed Apr 28, 2018
1 parent 40366a6 commit 289812b
Show file tree
Hide file tree
Showing 10 changed files with 2,880 additions and 15 deletions.
93 changes: 93 additions & 0 deletions src/lib/f_func.ts
Expand Up @@ -264,6 +264,8 @@ export interface Matrix {
lowerBand(value: number): Matrix;
packedUpper(value?: number): FortranArr;
packedLower(value?: number): FortranArr;
real(): Matrix;
imaginary(): Matrix;
toArr(): Complex[] | number[];
}

Expand Down Expand Up @@ -446,11 +448,102 @@ function mimicFMatrix(r: fpArray, i?: fpArray) {
}
}
return this;
},
real() {
let reC: fpArray;
reC = this.r.slice();
return mimicFMatrix(reC)(this.nrCols, this.colSize);
},
imaginary() {
let imC: fpArray;
if (!this.i) {
imC = new Float64Array(this.nrCols * this.colSize);
imC.fill(0);
}
else {
imC = this.i.slice(); //copy
}
return mimicFMatrix(imC)(this.nrCols, this.colSize);
}
});
}
}

export function real(data: number | number[] | Complex | Complex[]): number[] | number {

if (data === null) {
throw new Error('"null" inserted as argument');
}

if (data === undefined) {
throw new Error('"undefined inserted as argument');
}

if (typeof data === 'object' && 're' in data) {
return data.re;
}
if (typeof data === 'number') {
return data;
}
if (data instanceof Array) {
const copy: number[] = new Array<number>(data.length);
for (let i = 0; i < data.length; i++) {
const c = data[i];
if (typeof c === 'number') {
copy[i] = c;
continue;
}
if (c === undefined || c === null) {
copy[i] = c;
continue;
}
if ('re' in c) {
copy[i] = c['re'];
continue;
}
copy[i] = NaN;

}
return copy;
}
throw new Error('should not be here');
}

export function imaginary(data: number | number[] | Complex | Complex[]): number[] | number {

if (data === null) {
throw new Error('"null" inserted as argument');
}

if (data === undefined) {
throw new Error('"undefined inserted as argument');
}

if (typeof data === 'object' && 'im' in data) {
return data.im;
}

if (typeof data === 'number') {
return 0;
}

if (data instanceof Array) {
const copy: number[] = new Array<number>(data.length);
for (let i = 0; i <= data.length; i++) {
if (typeof data[i] === 'number') {
copy[i] = 0;
} else if ('im' in (data[i] as any)) {
copy[i] = data[i]['im'];
}
else {
copy[i] = NaN;
}
}
return copy;
}
throw new Error('should not be here');
}

export function fortranMatrixComplex32(...rest: (Complex | Complex[])[]):
(nrRows: number, nrCols: number, matrixType?: MatrixType, rowBase?: number, colBase?: number) => Matrix {

Expand Down
25 changes: 16 additions & 9 deletions src/lib/l3/single/sgemm.ts
@@ -1,4 +1,4 @@
import { errWrongArg, Matrix } from '../../f_func';
import { errWrongArg, lowerChar, Matrix } from '../../f_func';

/* -- Jacob Bogers, 03/2008, JS port, jkfbogers@gmail.cmom
*> -- Written on 8-February-1989.
Expand Down Expand Up @@ -26,8 +26,8 @@ export function sgemm(
ldc: number): void {

// faster then String.toLowerCase()
const trA: 'n' | 't' | 'c' = String.fromCharCode(transA.charCodeAt(0) | 0X20) as any;
const trB: 'n' | 't' | 'c' = String.fromCharCode(transB.charCodeAt(0) | 0X20) as any;
const trA = lowerChar(transA);
const trB = lowerChar(transB);

const notA = trA === 'n';
const notB = trB === 'n';
Expand All @@ -36,15 +36,15 @@ export function sgemm(
//ncolA is never used, I checked in the original F-code
// also checked online http://www.netlib.org/lapack/explore-html/db/dc9/group__single__blas__level3_gafe51bacb54592ff5de056acabd83c260.html#gafe51bacb54592ff5de056acabd83c260
// ncolA is not used
//const ncolA = notA ? k : n;
// const ncolA = notA ? k : n;

const nrowB = notB ? k : n;

let info = 0;
if (!notA && trA !== 'c' && trA !== 't') {
if (!'ntc'.includes(trA)) {
info = 1;
}
else if (!notB && trB !== 'c' && trB !== 't') {
else if (!'ntc'.includes(trB)) {
info = 2;
}
else if (m < 0) {
Expand All @@ -62,7 +62,7 @@ export function sgemm(
else if (ldb < max(1, nrowB)) {
info = 10;
}
else if (ldb < max(1, m)) {
else if (ldc < max(1, m)) {
info = 13;
}
// ok?
Expand All @@ -80,11 +80,18 @@ export function sgemm(

if (alpha === 0) {
if (beta === 0) {
c.r.fill(0); //fast
for (let j = 1; j <= n; j++) {
c.setCol(j, 1, m, 0);
}
}
else {
// I have to typecast it this way for TS compiler not to nag!!
(c.r as Float32Array).set((c.r as Float32Array).map(v => v * beta), 0);
for (let j = 1; j <= n; j++) {
const coorCJ = c.colOfEx(j);
for (let i = 1; i <= m; i++) {
c.r[coorCJ + i] *= beta;
}
}
}
return;
}
Expand Down
2 changes: 2 additions & 0 deletions test/generate-fixtures/sgemm/Makefile
@@ -0,0 +1,2 @@
test.exe : dgemm.f tests.f ../_shared/lsame.f ../_shared/xerbla.f
gfortran dgemm.f tests.f ../_shared/lsame.f ../_shared/xerbla.f -o test.exe

0 comments on commit 289812b

Please sign in to comment.