org.ejml.dense.block.InnerRankUpdate_FDRB Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ejml-fdense Show documentation
Show all versions of ejml-fdense Show documentation
A fast and easy to use dense and sparse matrix linear algebra library written in Java.
/*
* Copyright (c) 2009-2020, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.ejml.dense.block;
import javax.annotation.Generated;
import org.ejml.data.FMatrixRBlock;
import org.ejml.data.FSubmatrixD1;
//CONCURRENT_INLINE import static org.ejml.dense.block.InnerRankUpdate_FDRB.*;
//CONCURRENT_INLINE import org.ejml.concurrency.EjmlConcurrency;
/**
* Performs rank-n update operations on the inner blocks of a {@link FMatrixRBlock}
*
* It is assumed and not checked that the submatrices are aligned along the matrix's blocks.
*
* @author Peter Abeles
*/
@Generated("org.ejml.dense.block.InnerRankUpdate_DDRB")
public class InnerRankUpdate_FDRB {
/**
*
* Performs:
*
* A = A + α B TB
*
*
* @param blockLength Size of the block in the block matrix.
* @param alpha scaling factor for right hand side.
* @param A Block aligned submatrix.
* @param B Block aligned submatrix.
*/
public static void rankNUpdate( int blockLength, float alpha, FSubmatrixD1 A, FSubmatrixD1 B ) {
int heightB = B.row1 - B.row0;
if (heightB > blockLength)
throw new IllegalArgumentException("Height of B cannot be greater than the block length");
int N = B.col1 - B.col0;
if (A.col1 - A.col0 != N)
throw new IllegalArgumentException("A does not have the expected number of columns based on B's width");
if (A.row1 - A.row0 != N)
throw new IllegalArgumentException("A does not have the expected number of rows based on B's width");
//CONCURRENT_BELOW EjmlConcurrency.loopFor(B.col0,B.col1,blockLength,i->{
for (int i = B.col0; i < B.col1; i += blockLength) {
int indexB_i = B.row0*B.original.numCols + i*heightB;
int widthB_i = Math.min(blockLength, B.col1 - i);
int rowA = i - B.col0 + A.row0;
int heightA = Math.min(blockLength, A.row1 - rowA);
for (int j = B.col0; j < B.col1; j += blockLength) {
int widthB_j = Math.min(blockLength, B.col1 - j);
int indexA = rowA*A.original.numCols + (j - B.col0 + A.col0)*heightA;
int indexB_j = B.row0*B.original.numCols + j*heightB;
InnerMultiplication_FDRB.blockMultPlusTransA(alpha,
B.original.data, B.original.data, A.original.data,
indexB_i, indexB_j, indexA, heightB, widthB_i, widthB_j);
}
}
//CONCURRENT_ABOVE });
}
/**
*
* Rank N update function for a symmetric inner submatrix and only operates on the upper
* triangular portion of the submatrix.
*
* A = A - B TB
*
*/
public static void symmRankNMinus_U( int blockLength,
FSubmatrixD1 A, FSubmatrixD1 B ) {
int heightB = B.row1 - B.row0;
if (heightB > blockLength)
throw new IllegalArgumentException("Height of B cannot be greater than the block length");
int N = B.col1 - B.col0;
if (A.col1 - A.col0 != N)
throw new IllegalArgumentException("A does not have the expected number of columns based on B's width");
if (A.row1 - A.row0 != N)
throw new IllegalArgumentException("A does not have the expected number of rows based on B's width");
//CONCURRENT_BELOW EjmlConcurrency.loopFor(B.col0,B.col1,blockLength,i->{
for (int i = B.col0; i < B.col1; i += blockLength) {
int indexB_i = B.row0*B.original.numCols + i*heightB;
int widthB_i = Math.min(blockLength, B.col1 - i);
int rowA = i - B.col0 + A.row0;
int heightA = Math.min(blockLength, A.row1 - rowA);
for (int j = i; j < B.col1; j += blockLength) {
int widthB_j = Math.min(blockLength, B.col1 - j);
int indexA = rowA*A.original.numCols + (j - B.col0 + A.col0)*heightA;
int indexB_j = B.row0*B.original.numCols + j*heightB;
if (i == j) {
// only the upper portion of this block needs to be modified since it is along a diagonal
multTransABlockMinus_U(B.original.data, A.original.data,
indexB_i, indexB_j, indexA, heightB, widthB_i, widthB_j);
} else {
multTransABlockMinus(B.original.data, A.original.data,
indexB_i, indexB_j, indexA, heightB, widthB_i, widthB_j);
}
}
}
//CONCURRENT_ABOVE });
}
/**
*
* Rank N update function for a symmetric inner submatrix and only operates on the lower
* triangular portion of the submatrix.
*
* A = A - B*BT
*
*/
public static void symmRankNMinus_L( int blockLength,
FSubmatrixD1 A, FSubmatrixD1 B ) {
int widthB = B.col1 - B.col0;
if (widthB > blockLength)
throw new IllegalArgumentException("Width of B cannot be greater than the block length");
int N = B.row1 - B.row0;
if (A.col1 - A.col0 != N)
throw new IllegalArgumentException("A does not have the expected number of columns based on B's height");
if (A.row1 - A.row0 != N)
throw new IllegalArgumentException("A does not have the expected number of rows based on B's height");
//CONCURRENT_BELOW EjmlConcurrency.loopFor(B.row0,B.row1,blockLength,i->{
for (int i = B.row0; i < B.row1; i += blockLength) {
int heightB_i = Math.min(blockLength, B.row1 - i);
int indexB_i = i*B.original.numCols + heightB_i*B.col0;
int rowA = i - B.row0 + A.row0;
int heightA = Math.min(blockLength, A.row1 - rowA);
for (int j = B.row0; j <= i; j += blockLength) {
int widthB_j = Math.min(blockLength, B.row1 - j);
int indexA = rowA*A.original.numCols + (j - B.row0 + A.col0)*heightA;
int indexB_j = j*B.original.numCols + widthB_j*B.col0;
if (i == j) {
multTransBBlockMinus_L(B.original.data, A.original.data,
indexB_i, indexB_j, indexA, widthB, heightB_i, widthB_j);
} else {
multTransBBlockMinus(B.original.data, A.original.data,
indexB_i, indexB_j, indexA, widthB, heightB_i, widthB_j);
}
}
}
//CONCURRENT_ABOVE });
}
//CONCURRENT_OMIT_BEGIN
/**
*
* Performs the following operation on a block:
*
* c = c - aTa
*
*/
protected static void multTransABlockMinus( float[] dataA, float[] dataC,
int indexA, int indexB, int indexC,
final int heightA, final int widthA, final int widthC ) {
// for (int i = 0; i < widthA; i++) {
// for (int k = 0; k < heightA; k++) {
//
// float valA = dataA[k*widthA + i + indexA];
// for (int j = 0; j < widthC; j++) {
// dataC[i*widthC + j + indexC] -= valA*dataA[k*widthC + j + indexB];
// }
// }
// }
int rowB = indexB;
int endLoopK = rowB + heightA*widthC;
int startA = indexA;
//for( int k = 0; k < heightA; k++ ) {
for (; rowB != endLoopK; rowB += widthC, startA += widthA) {
int a = startA;
int c = indexC;
int endA = a + widthA;
int endB = rowB + widthC;
while (a != endA) {
float valA = dataA[a++];
int b = rowB;
while (b != endB) {
dataC[c++] -= valA*dataA[b++];
}
}
}
}
/**
*
* Performs the following operation on the upper triangular portion of a block:
*
* c = c - aTa
*
*/
protected static void multTransABlockMinus_U( float[] dataA, float[] dataC,
int indexA, int indexB, int indexC,
final int heightA, final int widthA, final int widthC ) {
// for (int i = 0; i < widthA; i++) {
// for (int k = 0; k < heightA; k++) {
//
// float valA = dataA[k*widthA + i + indexA];
// for (int j = i; j < widthC; j++) {
// dataC[i*widthC + j + indexC] -= valA*dataA[k*widthC + j + indexB];
// }
// }
// }
for (int i = 0; i < widthA; i++) {
for (int k = 0; k < heightA; k++) {
float valA = dataA[k*widthA + i + indexA];
int b = k*widthC + indexB + i;
int c = i*widthC + indexC + i;
int endC = (c - i) + widthC;
while (c != endC) {
// for( int j = i; j < widthC; j++ ) {
dataC[c++] -= valA*dataA[b++];
}
}
}
}
/**
*
* Performs the following operation on a block:
*
* c = c - a*aT
*
*/
protected static void multTransBBlockMinus( final float[] dataA, final float[] dataC,
final int indexA, final int indexB, final int indexC,
final int widthA, final int heightA, final int widthC ) {
// for (int i = 0; i < heightA; i++) {
// for (int j = 0; j < widthC; j++) {
// float sum = 0;
// for (int k = 0; k < widthA; k++) {
// sum += dataA[i*widthA + k + indexA]*dataA[j*widthA + k + indexB];
// }
// dataC[i*widthC + j + indexC] -= sum;
// }
// }
int rowA = indexA;
int c = indexC;
for (int i = 0; i < heightA; i++, rowA += widthA) {
final int endA = rowA + widthA;
int rowB = indexB;
final int endLoopJ = c + widthC;
// for( int j = 0; j < widthC; j++ ) {
while (c != endLoopJ) {
int a = rowA;
int b = rowB;
float sum = 0;
while (a != endA) {
sum += dataA[a++]*dataA[b++];
}
dataC[c++] -= sum;
rowB += widthA;
}
}
}
/**
*
* Performs the following operation on the lower triangular portion of a block:
*
* c = c - a*aT
*
*/
protected static void multTransBBlockMinus_L( float[] dataA, float[] dataC,
int indexA, int indexB, int indexC,
final int widthA, final int heightA, final int widthC ) {
// for (int i = 0; i < heightA; i++) {
// for (int j = 0; j <= i; j++) {
// float sum = 0;
// for (int k = 0; k < widthA; k++) {
// sum += dataA[i*widthA + k + indexA]*dataA[j*widthA + k + indexB];
// }
// dataC[i*widthC + j + indexC] -= sum;
// }
// }
for (int i = 0; i < heightA; i++) {
int rowA = i*widthA + indexA;
int endA = rowA + widthA;
int rowB = indexB;
int rowC = i*widthC + indexC;
for (int j = 0; j <= i; j++, rowB += widthA) {
float sum = 0;
int a = rowA;
int b = rowB;
while (a != endA) {
sum += dataA[a++]*dataA[b++];
}
dataC[rowC + j] -= sum;
}
}
}
//CONCURRENT_OMIT_END
}