diff --git a/CHANGELOG.md b/CHANGELOG.md
index ad6e80c..e5bdfb1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
 # Change log
 
+## 0.10.1
+ * Feature: Accept vectors in `matmul` and `fractionalMatmul` and adjust return types accordingly.
+ * Feature: `matscale`, `matadd`, `matsub` and fractional variants. [#41](https://github.com/xenharmonic-devs/xen-dev-utils/issues/41)
 ## 0.10.0
  * Feature: New array type `FractionalMonzo` with vector methods `fractionalMonzosEqual`, `fractionalAdd`, `fractionalSub`, `fractionalDot`, `fractionalScale` and `fractionalNorm`.
  * Feature: Unnormalized Gram-Schmidt orthogonalization `gram` and `fractionalGram`.
diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts
index 2f8a531..cc93e94 100644
--- a/src/__tests__/basis.spec.ts
+++ b/src/__tests__/basis.spec.ts
@@ -7,11 +7,17 @@ import {
   fractionalGram,
   fractionalInv,
   fractionalLenstraLenstraLovasz,
+  fractionalMatadd,
+  fractionalMatscale,
+  fractionalMatsub,
   fractionalMatmul,
   fractionalTranspose,
   gram,
   inv,
   lenstraLenstraLovasz,
+  matadd,
+  matscale,
+  matsub,
   matmul,
   minor,
   transpose,
@@ -368,6 +374,33 @@ describe('Matrix multiplication', () => {
     ]);
   });
 
+  it('multiplies a matrix with a vector', () => {
+    const A = [
+      [1, 2],
+      [3, 4],
+    ];
+    const v = [5, 6];
+    const u = matmul(A, v);
+    expect(u).toEqual([17, 39]);
+  });
+
+  it('multiplies a vector with a matrix', () => {
+    const A = [
+      [1, 2],
+      [3, 4],
+    ];
+    const v = [5, 6];
+    const u = matmul(v, A);
+    expect(u).toEqual([23, 34]);
+  });
+
+  it('multiplies a vector with a vector', () => {
+    const u = [1, 2];
+    const v = [5, 6];
+    const s = matmul(u, v);
+    expect(s).toEqual(17);
+  });
+
   it('multiplies two matrices (fractions)', () => {
     const A = [
       [1, 0, 1],
@@ -388,6 +421,33 @@ describe('Matrix multiplication', () => {
       ['17/3', '9', '6'],
     ]);
   });
+
+  it('multiplies a matrix with a vector (fractions)', () => {
+    const A = [
+      [1, 2],
+      [3, 0.5],
+    ];
+    const v = [5, '1/3'];
+    const u = fractionalMatmul(A, v);
+    expect(u.map(f => f.toFraction())).toEqual(['17/3', '91/6']);
+  });
+
+  it('multiplies a vector with a matrix (fractions)', () => {
+    const A = [
+      [1, 2],
+      [3, '1/3'],
+    ];
+    const v = [5, 0.5];
+    const u = fractionalMatmul(v, A);
+    expect(u.map(f => f.toFraction())).toEqual(['13/2', '61/6']);
+  });
+
+  it('multiplies a vector with a vector (fractions)', () => {
+    const u = [0.5, 2];
+    const v = [5, '5/7'];
+    const s = fractionalMatmul(u, v);
+    expect(s.toFraction()).toBe('55/14');
+  });
 });
 
 describe('Matrix inverse', () => {
@@ -654,3 +714,64 @@ describe('Transpose', () => {
     ]);
   });
 });
+
+describe('Auxiliary matrix methods', () => {
+  it('scales', () => {
+    const A = [
+      [1, 2],
+      [3, 4],
+    ];
+    expect(matscale(A, 2)).toEqual([
+      [2, 4],
+      [6, 8],
+    ]);
+    expect(
+      fractionalMatscale(A, 2).map(row => row.map(f => f.valueOf()))
+    ).toEqual([
+      [2, 4],
+      [6, 8],
+    ]);
+  });
+
+  it('adds', () => {
+    const A = [
+      [1, 2],
+      [3, 4],
+    ];
+    const B = [
+      [-1, 5],
+      [6, 7],
+    ];
+    expect(matadd(A, B)).toEqual([
+      [0, 7],
+      [9, 11],
+    ]);
+    expect(
+      fractionalMatadd(A, B).map(row => row.map(f => f.valueOf()))
+    ).toEqual([
+      [0, 7],
+      [9, 11],
+    ]);
+  });
+
+  it('subtracts', () => {
+    const A = [
+      [1, 2],
+      [3, 4],
+    ];
+    const B = [
+      [-1, 5],
+      [6, 7],
+    ];
+    expect(matsub(A, B)).toEqual([
+      [2, -3],
+      [-3, -3],
+    ]);
+    expect(
+      fractionalMatsub(A, B).map(row => row.map(f => f.valueOf()))
+    ).toEqual([
+      [2, -3],
+      [-3, -3],
+    ]);
+  });
+});
diff --git a/src/basis.ts b/src/basis.ts
index 2829189..db5ce09 100644
--- a/src/basis.ts
+++ b/src/basis.ts
@@ -2,6 +2,8 @@ import {Fraction, FractionValue} from './fraction';
 import {
   FractionalMonzo,
   ProtoFractionalMonzo,
+  add,
+  fractionalAdd,
   fractionalDot,
   fractionalScale,
   fractionalSub,
@@ -437,12 +439,35 @@ export function fractionalInv(matrix: ProtoFractionalMonzo[]) {
 }
 
 /**
- * Compute the matrix product of two arrays of arrays of numbers.
+ * Compute the matrix product of two matrices or vectors.
  * @param A The left operand.
  * @param B The right operand.
  * @returns The matrix product of the operands.
  */
-export function matmul(A: number[][], B: number[][]) {
+export function matmul(A: number[], B: number[]): number;
+export function matmul(A: number[], B: number[][]): number[];
+export function matmul(A: number[][], B: number[]): number[];
+export function matmul(A: number[][], B: number[][]): number[][];
+export function matmul(A: number[] | number[][], B: number[] | number[][]) {
+  let numVectors = 0;
+  if (!Array.isArray(A[0])) {
+    A = [A as number[]];
+    numVectors++;
+  }
+  if (!Array.isArray(B[0])) {
+    B = (B as number[]).map(c => [c]);
+    numVectors++;
+  }
+  const result = matmul_(A as number[][], B as number[][]);
+  if (numVectors === 1) {
+    return result.flat();
+  } else if (numVectors === 2) {
+    return result[0][0];
+  }
+  return result;
+}
+
+function matmul_(A: number[][], B: number[][]) {
   const height = A.length;
   let width = 0;
   for (const row of B) {
@@ -471,14 +496,55 @@ export function matmul(A: number[][], B: number[][]) {
 }
 
 /**
- * Compute the matrix product of two arrays of arrays of fractions.
+ * Compute the matrix product of two matrices or vectors.
  * @param A The left operand.
  * @param B The right operand.
  * @returns The matrix product of the operands.
  */
+export function fractionalMatmul(
+  A: ProtoFractionalMonzo,
+  B: ProtoFractionalMonzo
+): Fraction;
+export function fractionalMatmul(
+  A: ProtoFractionalMonzo,
+  B: ProtoFractionalMonzo[]
+): FractionalMonzo;
+export function fractionalMatmul(
+  A: ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo
+): FractionalMonzo;
 export function fractionalMatmul(
   A: ProtoFractionalMonzo[],
   B: ProtoFractionalMonzo[]
+): FractionalMonzo[];
+export function fractionalMatmul(
+  A: ProtoFractionalMonzo | ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo | ProtoFractionalMonzo[]
+) {
+  let numVectors = 0;
+  if (!Array.isArray(A[0])) {
+    A = [A as ProtoFractionalMonzo];
+    numVectors++;
+  }
+  if (!Array.isArray(B[0])) {
+    B = (B as ProtoFractionalMonzo).map(c => [c]);
+    numVectors++;
+  }
+  const result = fractionalMatmul_(
+    A as ProtoFractionalMonzo[],
+    B as ProtoFractionalMonzo[]
+  );
+  if (numVectors === 1) {
+    return result.flat();
+  } else if (numVectors === 2) {
+    return result[0][0];
+  }
+  return result;
+}
+
+export function fractionalMatmul_(
+  A: ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo[]
 ): FractionalMonzo[] {
   const height = A.length;
   let width = 0;
@@ -682,3 +748,92 @@ export function minor(matrix: any[][], i: number, j: number) {
   }
   return matrix;
 }
+
+/**
+ * Scale a matrix by a scalar.
+ * @param matrix The matrix to scale.
+ * @param amount The amount to scale by.
+ * @returns The scalar multiple.
+ */
+export function matscale(matrix: number[][], amount: number) {
+  return matrix.map(row => scale(row, amount));
+}
+
+/**
+ * Add two matrices.
+ * @param A The first matrix.
+ * @param B The second matrix.
+ * @returns The sum.
+ */
+export function matadd(A: number[][], B: number[][]) {
+  const result: number[][] = [];
+  const numRows = Math.max(A.length, B.length);
+  for (let i = 0; i < numRows; ++i) {
+    result.push(add(A[i] ?? [], B[i] ?? []));
+  }
+  return result;
+}
+
+/**
+ * Subtract two matrices.
+ * @param A The matrix to subtract from.
+ * @param B The matrix to subtract by.
+ * @returns The difference.
+ */
+export function matsub(A: number[][], B: number[][]) {
+  const result: number[][] = [];
+  const numRows = Math.max(A.length, B.length);
+  for (let i = 0; i < numRows; ++i) {
+    result.push(sub(A[i] ?? [], B[i] ?? []));
+  }
+  return result;
+}
+
+/**
+ * Scale a matrix by a scalar.
+ * @param matrix The matrix to scale.
+ * @param amount The amount to scale by.
+ * @returns The scalar multiple.
+ */
+export function fractionalMatscale(
+  matrix: ProtoFractionalMonzo[],
+  amount: FractionValue
+) {
+  return matrix.map(row => fractionalScale(row, amount));
+}
+
+/**
+ * Add two matrices.
+ * @param A The first matrix.
+ * @param B The second matrix.
+ * @returns The sum.
+ */
+export function fractionalMatadd(
+  A: ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo[]
+) {
+  const result: FractionalMonzo[] = [];
+  const numRows = Math.max(A.length, B.length);
+  for (let i = 0; i < numRows; ++i) {
+    result.push(fractionalAdd(A[i] ?? [], B[i] ?? []));
+  }
+  return result;
+}
+
+/**
+ * Subtract two matrices.
+ * @param A The matrix to subtract from.
+ * @param B The matrix to subtract by.
+ * @returns The difference.
+ */
+export function fractionalMatsub(
+  A: ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo[]
+) {
+  const result: FractionalMonzo[] = [];
+  const numRows = Math.max(A.length, B.length);
+  for (let i = 0; i < numRows; ++i) {
+    result.push(fractionalSub(A[i] ?? [], B[i] ?? []));
+  }
+  return result;
+}