testcinvunit.pas

maths lib with source

        Inc(I);
    end;
end;


(*************************************************************************
Matrix multiplication
*************************************************************************)
procedure InternalMatrixMatrixMultiply(const A : TComplex2DArray;
     AI1 : Integer;
     AI2 : Integer;
     AJ1 : Integer;
     AJ2 : Integer;
     TransA : Boolean;
     const B : TComplex2DArray;
     BI1 : Integer;
     BI2 : Integer;
     BJ1 : Integer;
     BJ2 : Integer;
     TransB : Boolean;
     var C : TComplex2DArray;
     CI1 : Integer;
     CI2 : Integer;
     CJ1 : Integer;
     CJ2 : Integer);
var
    ARows : Integer;
    ACols : Integer;
    BRows : Integer;
    BCols : Integer;
    CRows : Integer;
    CCols : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    L : Integer;
    R : Integer;
    V : Complex;
    WORK : TComplex1DArray;
    Beta : Complex;
    Alpha : Complex;
    i_ : Integer;
    i1_ : Integer;
begin
    
    //
    // Pre-setup
    //
    K := Max(AI2-AI1+1, AJ2-AJ1+1);
    K := Max(K, BI2-BI1+1);
    K := Max(K, BJ2-BJ1+1);
    SetLength(WORK, K+1);
    Beta := C_Complex(0);
    Alpha := C_Complex(1);
    
    //
    // Setup
    //
    if  not TransA then
    begin
        ARows := AI2-AI1+1;
        ACols := AJ2-AJ1+1;
    end
    else
    begin
        ARows := AJ2-AJ1+1;
        ACols := AI2-AI1+1;
    end;
    if  not TransB then
    begin
        BRows := BI2-BI1+1;
        BCols := BJ2-BJ1+1;
    end
    else
    begin
        BRows := BJ2-BJ1+1;
        BCols := BI2-BI1+1;
    end;
    Assert(ACols=BRows, 'MatrixMatrixMultiply: incorrect matrix sizes!');
    if (ARows<=0) or (ACols<=0) or (BRows<=0) or (BCols<=0) then
    begin
        Exit;
    end;
    CRows := ARows;
    CCols := BCols;
    
    //
    // Test WORK
    //
    I := Max(ARows, ACols);
    I := Max(BRows, I);
    I := Max(I, BCols);
    Work[1] := C_Complex(0);
    Work[I] := C_Complex(0);
    
    //
    // Prepare C
    //
    if C_EqualR(Beta,0) then
    begin
        I:=CI1;
        while I<=CI2 do
        begin
            J:=CJ1;
            while J<=CJ2 do
            begin
                C[I,J] := C_Complex(0);
                Inc(J);
            end;
            Inc(I);
        end;
    end
    else
    begin
        I:=CI1;
        while I<=CI2 do
        begin
            for i_ := CJ1 to CJ2 do
            begin
                C[I,i_] := C_Mul(Beta, C[I,i_]);
            end;
            Inc(I);
        end;
    end;
    
    //
    // A*B
    //
    if  not TransA and  not TransB then
    begin
        L:=AI1;
        while L<=AI2 do
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                V := C_Mul(Alpha,A[L,AJ1+R-BI1]);
                K := CI1+L-AI1;
                i1_ := (BJ1) - (CJ1);
                for i_ := CJ1 to CJ2 do
                begin
                    C[K,i_] := C_Add(C[K,i_], C_Mul(V, B[R,i_+i1_]));
                end;
                Inc(R);
            end;
            Inc(L);
        end;
        Exit;
    end;
    
    //
    // A*B'
    //
    if  not TransA and TransB then
    begin
        if ARows*ACols<BRows*BCols then
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                L:=AI1;
                while L<=AI2 do
                begin
                    i1_ := (BJ1)-(AJ1);
                    V := C_Complex(0.0);
                    for i_ := AJ1 to AJ2 do
                    begin
                        V := C_Add(V,C_Mul(A[L,i_],B[R,i_+i1_]));
                    end;
                    C[CI1+L-AI1,CJ1+R-BI1] := C_Add(C[CI1+L-AI1,CJ1+R-BI1],C_Mul(Alpha,V));
                    Inc(L);
                end;
                Inc(R);
            end;
            Exit;
        end
        else
        begin
            L:=AI1;
            while L<=AI2 do
            begin
                R:=BI1;
                while R<=BI2 do
                begin
                    i1_ := (BJ1)-(AJ1);
                    V := C_Complex(0.0);
                    for i_ := AJ1 to AJ2 do
                    begin
                        V := C_Add(V,C_Mul(A[L,i_],B[R,i_+i1_]));
                    end;
                    C[CI1+L-AI1,CJ1+R-BI1] := C_Add(C[CI1+L-AI1,CJ1+R-BI1],C_Mul(Alpha,V));
                    Inc(R);
                end;
                Inc(L);
            end;
            Exit;
        end;
    end;
    
    //
    // A'*B
    //
    if TransA and  not TransB then
    begin
        L:=AJ1;
        while L<=AJ2 do
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                V := C_Mul(Alpha,A[AI1+R-BI1,L]);
                K := CI1+L-AJ1;
                i1_ := (BJ1) - (CJ1);
                for i_ := CJ1 to CJ2 do
                begin
                    C[K,i_] := C_Add(C[K,i_], C_Mul(V, B[R,i_+i1_]));
                end;
                Inc(R);
            end;
            Inc(L);
        end;
        Exit;
    end;
    
    //
    // A'*B'
    //
    if TransA and TransB then
    begin
        if ARows*ACols<BRows*BCols then
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                I:=1;
                while I<=CRows do
                begin
                    WORK[I] := C_Complex(0.0);
                    Inc(I);
                end;
                L:=AI1;
                while L<=AI2 do
                begin
                    V := C_Mul(Alpha,B[R,BJ1+L-AI1]);
                    K := CJ1+R-BI1;
                    i1_ := (AJ1) - (1);
                    for i_ := 1 to CRows do
                    begin
                        WORK[i_] := C_Add(WORK[i_], C_Mul(V, A[L,i_+i1_]));
                    end;
                    Inc(L);
                end;
                i1_ := (1) - (CI1);
                for i_ := CI1 to CI2 do
                begin
                    C[i_,K] := C_Add(C[i_,K], WORK[i_+i1_]);
                end;
                Inc(R);
            end;
            Exit;
        end
        else
        begin
            L:=AJ1;
            while L<=AJ2 do
            begin
                K := AI2-AI1+1;
                i1_ := (AI1) - (1);
                for i_ := 1 to K do
                begin
                    WORK[i_] := A[i_+i1_,L];
                end;
                R:=BI1;
                while R<=BI2 do
                begin
                    i1_ := (BJ1)-(1);
                    V := C_Complex(0.0);
                    for i_ := 1 to K do
                    begin
                        V := C_Add(V,C_Mul(WORK[i_],B[R,i_+i1_]));
                    end;
                    C[CI1+L-AJ1,CJ1+R-BI1] := C_Add(C[CI1+L-AJ1,CJ1+R-BI1],C_Mul(Alpha,V));
                    Inc(R);
                end;
                Inc(L);
            end;
            Exit;
        end;
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testcinvunit_test_silent():Boolean;
begin
    Result := TestCInv(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testcinvunit_test():Boolean;
begin
    Result := TestCInv(False);
end;


end.