Announcement

Collapse
No announcement yet.

Need code to do matrix multiplication

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Need code to do matrix multiplication

    Hi All

    Is there some kinda PB code that can do matrix multiplication like the method use to multiply 2 matrices as given
    in https://www.mathsisfun.com/algebra/m...ltiplying.html






  • #2
    Click image for larger version

Name:	matrix-code.jpg?resize=800,450.jpg
Views:	109
Size:	87.0 KB
ID:	795265
    <b>George W. Bleck</b>
    <img src='http://www.blecktech.com/myemail.gif'>

    Comment


    • #3
      George not sure what is your graphics? what I need is some kinda of code to do the following row of array1 * column of array2 and so on

      Click image for larger version

Name:	Matrix M.PNG
Views:	108
Size:	45.0 KB
ID:	795268

      Comment


      • #4
        Check Help
        MAT statement
        Purpose
        To simplify Matrix Algebra calculations.

        Syntax
        MAT a1() = CON 'Set all elements of a1() to one
        MAT a1() = CON(expr) 'Set all elements of a1() to value of expr
        MAT a1() = IDN 'Establish a1() as an identity matrix
        MAT a1() = ZER 'Set all elements of a1() to zero
        MAT a1() = a2() + a3() 'Addition
        MAT a1() = a2() 'Assignment
        MAT a1() = INV(a2()) 'Inversion
        MAT a1() = (expr) * a2() 'Scalar Multiplication
        MAT a1() = a2() - a3() 'Subtraction
        MAT a1() = a2() * a3() 'Multiplication
        MAT a1() = TRN(a2()) 'Transposition

        Comment


        • #5
          take a look at the MAT function in PBHelp.
          Walt Decker

          Comment


          • #6
            Originally posted by Tim Lakinir View Post
            George not sure what is your graphics?
            It's a joke, It's a recurrent image in the movie "The Matrix"


            Comment


            • #7
              I've written and developed a detailed matrix-vector-scalar library I'd be willing to post. It applies a specific vector and matrix type and works by pointer. The types are as follows:

              Code:
              ' UDT for VEC (vector) structure
              ' 0 to 2 are the ijk elements, with the 3rd element as the magnitude
              ' Modified 20160624, performed the conversion to true three element vector
              %i_VEC_x = 00
              %i_VEC_y = 01
              %i_VEC_z = 02
              
              ' Supports unit vectors
              %e_TYPE_UNIT_VECTOR_i = 01
              %e_TYPE_UNIT_VECTOR_j = 02
              %e_TYPE_UNIT_VECTOR_k = 03
              
              TYPE UDT_VEC
              
              v(00 TO 02) AS DOUBLE
              
              END TYPE

              Comment


              • #8
                you can gather type of functions in each lib:

                Code:
                #INCLUDE ONCE "SUB Matrix_ADD.bas"
                #INCLUDE ONCE "SUB Matrix_SUB.bas"
                #INCLUDE ONCE "SUB Matrix_MLT.bas"
                #INCLUDE ONCE "SUB Matrix_TRN.bas"
                #INCLUDE ONCE "SUB Matrix_IDN.bas"
                #INCLUDE ONCE "SUB Matrix_CPY.bas"
                #INCLUDE ONCE "SUB Matrix_MSGBOX.bas"
                #INCLUDE ONCE "SUB Matrix_R.bas"
                #INCLUDE ONCE "SUB Matrix_T_EAR.bas"
                #INCLUDE ONCE "SUB Matrix_T_RQF_xyz.bas"
                #INCLUDE ONCE "SUB Matrix_T_xyz_RQF.bas"
                #INCLUDE ONCE "SUB Matrix_T_RQF_UVW_MOD.bas"
                #INCLUDE ONCE "SUB Matrix_T_xyz_UVW_MOD.bas"
                #INCLUDE ONCE "SUB Matrix_T_xyz_uvw.bas"
                #INCLUDE ONCE "SUB Matrix_T_xyz_rst.bas"
                #INCLUDE ONCE "SUB Matrix_T_RQF_zsh.bas"
                #INCLUDE ONCE "SUB Matrix_U_xyz.bas"
                #INCLUDE ONCE "SUB Matrix_T_SGP_xyz.bas"
                
                #INCLUDE ONCE "SUB Vector_ADD.bas"
                #INCLUDE ONCE "SUB Vector_INT.bas"
                #INCLUDE ONCE "SUB Vector_DER.bas"
                #INCLUDE ONCE "SUB Vector_SUB.bas"
                #INCLUDE ONCE "SUB Vector_MLT.bas"
                #INCLUDE ONCE "SUB Vector_AVG.bas"
                #INCLUDE ONCE "SUB Vector_ABS.bas"
                #INCLUDE ONCE "SUB Vector_CRS.bas"
                #INCLUDE ONCE "SUB Vector_CPY.bas"
                #INCLUDE ONCE "SUB Vector_NORM.bas"
                #INCLUDE ONCE "SUB Vector_MSGBOX.bas"
                #INCLUDE ONCE "SUB Vector_INPUT.bas"
                #INCLUDE ONCE "SUB Vector_INTP.bas"
                #INCLUDE ONCE "SUB Vector_TAB_INTP.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_c_RQF.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_r_xyz.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_r_RQF.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_a_RQF_a_xyz.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_v_RQF_v_LLA.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_a_xyz_a_RQF.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_v_xyz_c_RQF_t.bas"
                #INCLUDE ONCE "SUB Vector_v_xyz_v_RQF_VVH"
                #INCLUDE ONCE "SUB Vector_r_xyz_c_RAD.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_r_xyz_BCF.bas"
                #INCLUDE ONCE "SUB Vector_c_RAD_r_xyz.bas"
                #INCLUDE ONCE "SUB Vector_c_RAD_c_RQF.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_c_RAD.bas"
                #INCLUDE ONCE "SUB Vector_c_RAD_c_RAD_2p.bas"
                #INCLUDE ONCE "SUB Vector_c_LLA_2p_c_LLA.bas"
                #INCLUDE ONCE "SUB Vector_c_LLA_v_LLA_c_LLA_t.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_c_LLA_t_v_RQF.bas"
                #INCLUDE ONCE "SUB Vector_c_LLA_c_LLA_t_v_LLA.bas"
                #INCLUDE ONCE "SUB Vector_c_LLA_c_RQF.bas"
                #INCLUDE ONCE "SUB Vector_c_LLA_r_xyz.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_c_LLA.bas"
                #INCLUDE ONCE "SUB Vector_c_RQF_c_LLA.bas"
                #INCLUDE ONCE "SUB Vector_ARW.bas"
                #INCLUDE ONCE "SUB Vector_UNIT.bas"
                #INCLUDE ONCE "SUB Vector_CONS.bas"
                #INCLUDE ONCE "SUB Vector_SET.bas"
                #INCLUDE ONCE "SUB Vector_r_xyz_v_xyz_U_KOE.bas"
                #INCLUDE ONCE "SUB Vector_U_KOE_r_xyz_v_xyz.bas"
                #INCLUDE ONCE "SUB Vector_xyz_EQT_XYZ_ECL.bas"
                #INCLUDE ONCE "SUB Vector_XYZ_ECL_xyz_EQT.bas"
                #INCLUDE ONCE "SUB Vector_RAD_EQT_RAD_ECL.bas"
                #INCLUDE ONCE "SUB Vector_RAD_ECL_RAD_EQT.bas"
                
                #INCLUDE ONCE "FCN Scalar_DET.bas"
                #INCLUDE ONCE "FCN Scalar_NORM.bas"
                #INCLUDE ONCE "FCN Scalar_DOT.bas"
                #INCLUDE ONCE "FCN Scalar_ARC.bas"
                #INCLUDE ONCE "FCN Scalar_ANG.bas"
                #INCLUDE ONCE "FCN Scalar_AREA.bas"
                #INCLUDE ONCE "FCN Scalar_c_LLA_c_LLA_t_v_LLA.bas"

                Comment


                • #9
                  As the matrix evolved, I could start to use it to build other functions. Here is an example of a variable Euler angle sequenced rotation. Basically you can define any permutation, like PYR or RPY or YPR etc.

                  Code:
                  SUB Matrix_T_EAR (w_EAR_ijk AS LONG , _
                  BYVAL pa_EAR AS UDT_VEC POINTER , _
                  BYVAL pM_DCM AS DOUBLE POINTER )
                  
                  ' Routine performs a sequenced Euler angle rotation (EAR)
                  
                  DIM w_EAR_i AS LONG
                  DIM w_EAR_j AS LONG
                  DIM w_EAR_k AS LONG
                  
                  DIM pR_i AS DWORD
                  DIM pR_j AS DWORD
                  DIM pR_k AS DWORD
                  
                  DIM L_INIT AS STATIC LONG
                  
                  IF NOT L_INIT THEN
                  L_INIT = %T
                  REDIM M_i (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                  REDIM M_j (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                  REDIM M_ji(00 TO 02, 00 TO 02) AS STATIC DOUBLE
                  REDIM M_k (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                  END IF
                  
                  
                  IF pa_EAR = 00 OR _
                  pM_DCM = 00 THEN
                  
                  EXIT SUB
                  
                  ELSE
                  
                  ' Test for default condition which is the standard 321 sequence
                  ' Note this is done to speed the most often used 321 sequence
                  
                  IF w_EAR_ijk = 321 THEN
                  
                  ' Creates a DCM by Euler angle rotation 321
                  
                  @pM_DCM[00 OF 02, 00 OF 02] = COS(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                  @pM_DCM[00 OF 02, 01 OF 02] = COS(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                  @pM_DCM[00 OF 02, 02 OF 02] = - SIN(@pa_EAR.v(01))
                  
                  @pM_DCM[01 OF 02, 00 OF 02] = - COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(00)) + SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                  @pM_DCM[01 OF 02, 01 OF 02] = COS(@pa_EAR.v(02)) * COS(@pa_EAR.v(00)) + SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                  @pM_DCM[01 OF 02, 02 OF 02] = SIN(@pa_EAR.v(02)) * COS(@pa_EAR.v(01))
                  
                  @pM_DCM[02 OF 02, 00 OF 02] = SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(00)) + COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                  @pM_DCM[02 OF 02, 01 OF 02] = - SIN(@pa_EAR.v(02)) * COS(@pa_EAR.v(00)) + COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                  @pM_DCM[02 OF 02, 02 OF 02] = COS(@pa_EAR.v(02)) * COS(@pa_EAR.v(01))
                  
                  
                  ELSE
                  
                  ' Parse out the EAR sequence ijk into individual components
                  ' Note this will handle any combination including zeros which are trapped below
                  w_EAR_i = FIX( w_EAR_ijk / 100.0#)
                  w_EAR_j = FIX((w_EAR_ijk - 100 * w_EAR_i) / 10.0#)
                  w_EAR_k = FIX((w_EAR_ijk - 100 * w_EAR_i - 10 * w_EAR_j) / 1.0#)
                  
                  #IF 00
                  MSGBOX "w_EAR i,j,k = " + STR$(w_EAR_i) + " " + _
                  STR$(w_EAR_j) + " " + _
                  STR$(w_EAR_k) , %MB_ICONINFORMATION OR %MB_TASKMODAL , $S_MSGBOX_HDR
                  #ENDIF
                  
                  ' Choose the codepointer for the specific rotation based on the EAR ijk indices
                  pR_i = CHOOSE(w_EAR_i, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                  pR_j = CHOOSE(w_EAR_j, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                  pR_k = CHOOSE(w_EAR_k, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                  
                  ' Test the function pointer and call the respective rotation matrix by code pointer
                  ' Note the null rotation pointer is trapped by applying the identity matrix
                  IF pR_i = 00 THEN
                  Matrix_IDN VARPTR(M_i(00))
                  ELSE
                  CALL DWORD pR_i USING Matrix_R_n_PATTERN(@pa_EAR.v(00), VARPTR(M_i(00)))
                  END IF
                  
                  IF pR_j = 00 THEN
                  Matrix_IDN VARPTR(M_j(00))
                  ELSE
                  CALL DWORD pR_j USING Matrix_R_n_PATTERN(@pa_EAR.v(01), VARPTR(M_j(00)))
                  END IF
                  
                  IF pR_k = 00 THEN
                  Matrix_IDN VARPTR(M_k(00))
                  ELSE
                  CALL DWORD pR_k USING Matrix_R_n_PATTERN(@pa_EAR.v(02), VARPTR(M_k(00)))
                  END IF
                  
                  ' Create Mji = Mj * Mi
                  Matrix_MLT VARPTR(M_j (00)) , _
                  VARPTR(M_i (00)) , _
                  VARPTR(M_ji(00))
                  
                  ' Create M_DCM = Mk * Mji = Mk * Mj * Mi
                  Matrix_MLT VARPTR(M_k (00)) , _
                  VARPTR(M_ji(00)) , _
                  pM_DCM
                  
                  END IF
                  
                  END IF
                  
                  END SUB
                  
                  #IF %F ' Used for testing
                  
                  REDIM M_EAR(00 TO 02, 00 TO 02) AS DOUBLE
                  REDIM a_EAR(00 TO 02) AS DOUBLE
                  
                  ' Create a DCM Euler matrix based on a ijk rotation
                  
                  a_EAR(00) = 10.0# / 180.0# * Number_PI
                  a_EAR(01) = 20.0# / 180.0# * Number_PI
                  a_EAR(02) = 30.0# / 180.0# * Number_PI
                  
                  Matrix_R_DCM_EUL 123 , _
                  VARPTR(a_EAR(00)) , _
                  VARPTR(M_EAR(00))
                  
                  
                  ? Matrix_GET_S(VARPTR(M_EAR(00)))
                  
                  #ENDIF

                  Comment


                  • #10
                    If all you want is multiplication, then simply extend the following Euclidean 3x3 model

                    Code:
                    SUB Matrix_MLT ( _
                    BYVAL pA_MAT AS DOUBLE POINTER , _
                    BYVAL pB_MAT AS DOUBLE POINTER , _
                    BYVAL pC_MAT AS DOUBLE POINTER )
                    
                    ' Routine performs a left right matrix multiplication [C] = [A] [B]
                    
                    DIM i_ROW AS LONG
                    DIM j_COL AS LONG
                    
                    
                    IF pA_MAT = 00 OR _
                    pB_MAT = 00 OR _
                    pC_MAT = 00 THEN
                    
                    EXIT SUB
                    
                    ELSE
                    
                    FOR i_ROW = 00 TO 02
                    
                    FOR j_COL = 00 TO 02
                    
                    @pC_MAT[i_ROW OF 02, j_COL OF 02] = @pA_MAT[i_ROW OF 02, 00 OF 02] * @pB_MAT[00 OF 02, j_COL OF 02] + _
                    @pA_MAT[i_ROW OF 02, 01 OF 02] * @pB_MAT[01 OF 02, j_COL OF 02] + _
                    @pA_MAT[i_ROW OF 02, 02 OF 02] * @pB_MAT[02 OF 02, j_COL OF 02]
                    
                    NEXT j_COL
                    
                    NEXT i_ROW
                    
                    END IF
                    
                    END SUB

                    Comment


                    • #11
                      I extended this for variable row-column formats:

                      Code:
                      #COMPILE EXE
                      #DIM ALL
                      
                      FUNCTION PBMAIN () AS LONG
                      
                      REDIM A(00 TO 02 - 01, 00 TO 04 - 01) AS DOUBLE
                      REDIM B(00 TO 04 - 01, 00 TO 02 - 01) AS DOUBLE
                      REDIM C(00 TO 02 - 01, 00 TO 02 - 01) AS DOUBLE
                      
                      A(1 - 1, 1 - 1) = 1.1# : A(1 - 1, 2 - 1) = 1.2# : A(1 - 1, 3 - 1) = 1.3# : A(1 - 1, 4 - 1) = 1.4#
                      A(2 - 1, 1 - 1) = 2.1# : A(2 - 1, 2 - 1) = 2.2# : A(2 - 1, 3 - 1) = 2.3# : A(2 - 1, 4 - 1) = 2.4#
                      
                      B(1 - 1, 1 - 1) = 1.1# : B(1 - 1, 2 - 1) = 1.2#
                      B(2 - 1, 1 - 1) = 2.1# : B(2 - 1, 2 - 1) = 2.2#
                      B(3 - 1, 1 - 1) = 3.1# : B(3 - 1, 2 - 1) = 3.2#
                      B(4 - 1, 1 - 1) = 4.1# : B(4 - 1, 2 - 1) = 4.2#
                      
                      Matrix_MLT_VAR 02 , _
                      04 , _
                      04 , _
                      02 , _
                      VARPTR(A(00)) , _
                      VARPTR(B(00)) , _
                      VARPTR(C(00))
                      
                      ? "A" + $CRLF + _
                      Matrix_GET_S_VAR (02 , _
                      04 , _
                      VARPTR(A(00)))
                      
                      ? "B" + $CRLF + _
                      Matrix_GET_S_VAR (04 , _
                      02 , _
                      VARPTR(B(00)))
                      
                      ? "AxB" + $CRLF + _
                      Matrix_GET_S_VAR (02 , _
                      02 , _
                      VARPTR(C(00)))
                      
                      
                      END FUNCTION
                      
                      
                      SUB Matrix_MLT_VAR ( _
                      N_ROW_A AS LONG , _
                      N_COL_A AS LONG , _
                      N_ROW_B AS LONG , _
                      N_COL_B AS LONG , _
                      BYVAL pA_MAT AS DOUBLE POINTER , _
                      BYVAL pB_MAT AS DOUBLE POINTER , _
                      BYVAL pC_MAT AS DOUBLE POINTER )
                      
                      ' Routine performs a left right matrix multiplication [C] = [A] [B]
                      ' Note this is a variable row-column variant
                      
                      DIM i_ROW AS LONG
                      DIM j_COL AS LONG
                      
                      DIM k_INS AS LONG
                      DIM N_INS AS LONG
                      DIM N_OUT AS LONG
                      
                      
                      IF N_ROW_A = 00 OR _
                      N_COL_A = 00 OR _
                      N_ROW_B = 00 OR _
                      N_COL_B = 00 OR _
                      pA_MAT = 00 OR _
                      pB_MAT = 00 OR _
                      pC_MAT = 00 THEN
                      
                      EXIT SUB
                      
                      ELSEIF NOT N_COL_A = N_ROW_B THEN
                      
                      ' Trap improper row column form
                      
                      EXIT SUB
                      
                      ELSE
                      
                      N_INS = (N_COL_A + N_ROW_B) \ 02
                      N_OUT = (N_ROW_A + N_COL_B) \ 02
                      
                      FOR i_ROW = 01 TO N_OUT
                      
                      FOR j_COL = 01 TO N_OUT
                      
                      FOR k_INS = 01 TO N_INS
                      
                      @pC_MAT[i_ROW - 01 OF N_OUT - 01, j_COL - 01 OF N_OUT - 01] += @pA_MAT[i_ROW - 01 OF N_ROW_A - 01, k_INS - 01 OF N_COL_A - 01] * @pB_MAT[k_INS - 01 OF N_ROW_B - 01, j_COL - 01 OF N_COL_B - 01]
                      
                      NEXT k_INS
                      
                      NEXT j_COL
                      
                      NEXT i_ROW
                      
                      END IF
                      
                      END SUB
                      
                      
                      FUNCTION Matrix_GET_S_VAR ( _
                      N_ROW_A AS LONG , _
                      N_COL_A AS LONG , _
                      BYVAL pA_MAT AS DOUBLE POINTER ) AS STRING
                      
                      ' Routine returns the matrix as a string
                      ' Note this applies to a variable row column matrix
                      
                      DIM i_ROW AS LONG
                      DIM j_COL AS LONG
                      
                      DIM S_MAT AS STRING
                      
                      IF N_ROW_A = 00 OR _
                      N_COL_A = 00 OR _
                      pA_MAT = 00 THEN
                      
                      FUNCTION = "ERROR"
                      
                      ELSE
                      
                      FOR i_ROW = 01 TO N_ROW_A
                      
                      FOR j_COL = 01 TO N_COL_A
                      
                      S_MAT += STR$(@pA_MAT[i_ROW - 01 OF N_ROW_A - 01, j_COL - 01 OF N_COL_A - 01], 10) + IIF$(j_COL = N_COL_A, $CRLF, " ")
                      
                      NEXT j_COL
                      
                      NEXT i_ROW
                      
                      FUNCTION = S_MAT
                      
                      END IF
                      
                      END FUNCTION

                      Comment


                      • #12
                        Thanks Sir Dean, your code works very well
                        I had modified it a bit so that I can understand it better and I use the example given in https://www.mathsisfun.com/algebra/m...ltiplying.html
                        to test it and the answer is validated

                        Here's my rendition of your code

                        Code:
                        ' Matrix multiplication.bas
                        
                          ' Thanks to Dean Schrage
                        
                        #COMPILE EXE
                        #DIM ALL
                        #INCLUDE "win32api.inc"
                        
                        
                        '==============================
                        FUNCTION PBMAIN () AS LONG
                          LOCAL NumArow , NumAcol , NumBrow , NumBcol  AS LONG
                          LOCAL NumCrow , NumCcol AS LONG
                        
                          GOTO Exam2
                        
                         ' Example #1 -- original example by Dean
                           NumArow = 2
                           NumAcol = 4
                           NumBrow = 4
                           NumBcol = 2
                        
                           NumCrow= 2
                           NumCcol= 2
                        
                          REDIM A(1 TO 2 ,1 TO 4) AS DOUBLE
                          REDIM B(1 TO 4, 1 TO 2) AS DOUBLE
                          REDIM C(1 TO 2, 1 TO 2) AS DOUBLE
                        
                            A(1,1) = 1.1#
                            A(1,2) = 1.2#
                            A(1,3) = 1.3#
                            A(1,4) = 1.4#
                        
                            A(2,1) = 2.1#
                            A(2,2) = 2.2#
                            A(2,3) = 2.3#
                            A(2,4) = 2.4#
                        
                            B(1,1) = 1.1#
                            B(1,2) = 1.2#
                            B(2,1) = 2.1#
                            B(2,2) = 2.2#
                        
                        
                            B(3,1) = 3.1#
                            B(3,2) = 3.2#
                            B(4,1) = 4.1#
                            B(4,2) = 4.2#
                        
                        
                        
                        
                           Exam2 :
                          ' Example 2 from  https://www.mathsisfun.com/algebra/matrix-multiplying.html
                        
                           NumArow = 2
                           NumAcol = 3
                        
                           NumBrow = 3
                           NumBcol = 2
                        
                           NumCrow= 2
                           NumCcol= 2
                        
                           REDIM A(1 TO 2 ,1 TO 3) AS DOUBLE
                           REDIM B(1 TO 3, 1 TO 2) AS DOUBLE
                           REDIM C(1 TO 2, 1 TO 2) AS DOUBLE
                        
                            A(1,1) = 1#
                            A(1,2) = 2#
                            A(1,3) = 3#
                        
                            A(2,1) = 4#
                            A(2,2) = 5#
                            A(2,3) = 6#
                        
                            B(1,1) = 7#
                            B(1,2) = 8#
                            B(2,1) = 9#
                            B(2,2) = 10#
                        
                            B(3,1) = 11#
                            B(3,2) = 12#
                        
                        
                        
                        '  Do the matrix multiplication
                           Matrix_MLT_VAR( NumArow , NumAcol , NumBrow , NumBcol ,_
                                      VARPTR(A(1)) , VARPTR(B(1)) , VARPTR(C(1)))
                        
                        
                         ' resultant matrix strings
                           LOCAL AryA , AryB , AryC AS STRING
                        
                           AryA  = "[A] elements :" + $CRLF + _
                                Matrix_GET_S_VAR(NumArow , NumAcol ,VARPTR(A(1)))
                        
                           AryB = "[B] elements :" + $CRLF + _
                                Matrix_GET_S_VAR(NumBrow , NumBcol , VARPTR(B(1)))
                        
                           AryC = "results of [A] * [B]  elements :" + $CRLF + _
                              Matrix_GET_S_VAR(NumCrow , NumCcol , VARPTR(C(1)))
                        
                        
                         '  output to a text file
                            OPEN "Results.txt" FOR OUTPUT AS #3
                            PRINT #3, AryA
                            PRINT #3,
                            PRINT #3,
                        
                            PRINT #3, AryB
                            PRINT #3,
                            PRINT #3,
                        
                            PRINT #3, AryC
                            PRINT #3,
                            PRINT #3,
                        
                            CLOSE #3
                        
                        
                            ' Display output using notepad
                              LOCAL  lRes AS LONG
                              lRes = ShellExecute(0, "", "Notepad.exe", "Results.txt", "", %SW_SHOWNORMAL)
                        
                        END FUNCTION
                        
                        
                        
                        
                        '=========================
                        ' Matrix multiplication
                        ' Routine performs a left right matrix multiplication [C] = [A] * [B]
                        ' Note this is a variable row-column variant
                        SUB Matrix_MLT_VAR ( N_ARow AS LONG , N_ACol AS LONG , _
                              N_BRow AS LONG , N_BCol AS LONG , _
                              BYVAL pA_MAT AS DOUBLE POINTER , _
                              BYVAL pB_MAT AS DOUBLE POINTER , _
                              BYVAL pC_MAT AS DOUBLE POINTER )
                        
                        
                        
                        LOCAL i_Rw AS LONG
                        LOCAL j_Cl AS LONG
                        
                        LOCAL k_INS AS LONG
                        LOCAL N_INS AS LONG
                        LOCAL N_OUT AS LONG
                        
                        
                          IF N_ARow = 0 OR  N_ACol = 0 OR  N_BRow = 0 OR _
                             N_BCol = 0 OR  pA_MAT = 0 OR  pB_MAT = 0 OR _
                             pC_MAT = 0 THEN
                                 EXIT SUB
                          END IF
                        
                        
                           IF NOT N_ACol = N_BRow THEN
                        '     Trap improper row column form
                              EXIT SUB
                           END IF
                        
                          N_INS = (N_ACol + N_BRow) \ 2
                          N_OUT = (N_ARow + N_BCol) \ 2
                        
                          FOR i_Rw = 1 TO N_OUT
                              FOR j_Cl = 1 TO N_OUT
                        
                              FOR k_INS = 1 TO N_INS
                                  @pC_MAT[i_Rw - 1 OF N_OUT - 1, j_Cl - 1 OF N_OUT - 1] += _
                                  @pA_MAT[i_Rw - 1 OF N_ARow - 1, k_INS - 1 OF N_ACol - 1] * _
                                  @pB_MAT[k_INS - 1 OF N_BRow - 1, j_Cl - 1 OF N_BCol - 1]
                              NEXT k_INS
                        
                              NEXT j_Cl
                          NEXT i_Rw
                        
                        
                        END SUB
                        
                        
                        
                        '====================================
                        ' Routine returns the matrix as a string
                        ' Note this applies to a variable row column matrix
                        FUNCTION Matrix_GET_S_VAR ( N_Row AS LONG , N_Col AS LONG , _
                                      BYVAL ptr_Mat AS DOUBLE POINTER ) AS STRING
                        
                        
                        LOCAL i_Row AS LONG
                        LOCAL j_Col AS LONG
                        
                        LOCAL Str_Mat AS STRING
                        
                          IF N_Row = 0 OR N_Col = 0 OR ptr_Mat = 0 THEN
                              FUNCTION = "ERROR"
                          END IF
                        
                        
                          FOR i_Row = 1 TO N_Row
                            FOR j_Col = 1 TO N_Col
                                ' Matrix string
                                  Str_Mat += STR$(@ptr_Mat[i_Row - 1 OF N_Row - 1, j_Col - 1 OF N_Col - 1], 10) + _
                                         IIF$(j_Col = N_Col, $CRLF, " ")
                            NEXT j_Col
                          NEXT i_Row
                        
                          FUNCTION = Str_Mat
                        
                        
                        END FUNCTION

                        Comment


                        • #13
                          An other way Using PB MAT

                          Both example 1 and 2 of your link.

                          '
                          Code:
                          ' Matrix multiplication.bas
                          #COMPILE EXE
                          #DIM ALL
                          '==============================
                          FUNCTION PBMAIN () AS LONG
                          LOCAL AryA , AryB , AryC AS STRING
                          LOCAL S AS DOUBLE
                             REDIM B(1 TO 2, 1 TO 2) AS DOUBLE
                             REDIM C(1 TO 2, 1 TO 2) AS DOUBLE
                          
                          ' Example 1 from  https://www.mathsisfun.com/algebra/matrix-multiplying.html
                              S = 2
                          
                              B(1,1) = 4#
                              B(1,2) = 0#
                              B(2,1) = 1#
                              B(2,2) = -9#
                          
                          MAT C()= (S) * B()                        ' PB Matrix Scalar Multiplication
                          
                           ' resultant matrix strings
                             AryA  = "Single Number :" + $CRLF + _
                                  STR$(S) + $CRLF
                          
                             AryB = "[B] elements :" + $CRLF + _
                                  Matrix_GET_S_VAR(B())
                          
                             AryC = "results of S * [B]  elements :" + $CRLF + _
                                Matrix_GET_S_VAR(C())
                          
                          MSGBOX "Example 1" + $CRLF + AryA + $CRLF + AryB + $CRLF + AryC
                          
                          
                            ' Example 2 from  https://www.mathsisfun.com/algebra/matrix-multiplying.html
                          
                             REDIM A(1 TO 2 ,1 TO 3) AS DOUBLE
                             REDIM B(1 TO 3, 1 TO 2) AS DOUBLE
                             REDIM C(1 TO 2, 1 TO 2) AS DOUBLE
                          
                              A(1,1) = 1#
                              A(1,2) = 2#
                              A(1,3) = 3#
                          
                              A(2,1) = 4#
                              A(2,2) = 5#
                              A(2,3) = 6#
                          
                              B(1,1) = 7#
                              B(1,2) = 8#
                              B(2,1) = 9#
                              B(2,2) = 10#
                          
                              B(3,1) = 11#
                              B(3,2) = 12#
                          
                             MAT c()= a() * b()                 ' PB Matrix Multiply
                          
                           ' resultant matrix strings
                             AryA  = "[A] elements :" + $CRLF + _
                                  Matrix_GET_S_VAR(A())
                          
                             AryB = "[B] elements :" + $CRLF + _
                                  Matrix_GET_S_VAR(B())
                          
                             AryC = "results of [A] * [B]  elements :" + $CRLF + _
                                Matrix_GET_S_VAR(C())
                          
                          MSGBOX "Example 2" + $CRLF + AryA + $CRLF + AryB + $CRLF + AryC
                          END FUNCTION
                          
                          
                          '====================================
                          ' Routine returns the matrix as a string
                          ' Note this applies to a variable row column matrix
                          FUNCTION Matrix_GET_S_VAR (ARR() AS DOUBLE) AS STRING
                          
                          LOCAL i_Row, N_Row AS LONG
                          LOCAL j_Col, N_Col AS LONG
                          LOCAL Str_Mat AS STRING
                          LOCAL ptr_Mat AS DOUBLE POINTER
                          
                          ptr_Mat = VARPTR(ARR(1))
                          N_Row   = UBOUND(ARR(1))
                          N_Col   = UBOUND(ARR(2))
                          
                            FOR i_Row = 1 TO N_Row
                              FOR j_Col = 1 TO N_Col
                                  ' Matrix string
                                    Str_Mat += STR$(@ptr_Mat[i_Row - 1 OF N_Row - 1, j_Col - 1 OF N_Col - 1], 10) + _
                                           IIF$(j_Col = N_Col, $CRLF, " ")
                              NEXT j_Col
                            NEXT i_Row
                          
                            FUNCTION = Str_Mat
                          END FUNCTION
                          '

                          Comment


                          • #14
                            Well, nice to see it worked for you.

                            As for 1's based arrays, that used to be my standard but now I do a 0's based index since this is easier to derive pointers and pass with other languages. To each his own but if you want to use the library that I have you need a 0's based approach.

                            As for B MAT functions, you'll run aground eventually with this approach. I think it is much better to develop a generic library. Once I did this, I've ported it to C and Python. In fact using ctypes in Python (https://docs.python.org/3/library/ctypes.html), the library is just as fast as Numpy. So once you get a library, you are good to go. Here's a cross product and Vector norm both of which are used in OpenGL normal and lighting calls.

                            Code:
                            ' Process the orbit plane SAT CTR angular momentum vector in xyz ECI coordinates, note eh() is a normalized unit vector after the cross product
                            Vector_CRS VARPTR(r) , _
                            VARPTR(v) , _
                            VARPTR(eh)
                            
                            Scalar_NORM VARPTR(r) TO r_MAG
                            Scalar_NORM VARPTR(v) TO v_MAG
                            Code:
                            SUB Vector_CRS ( _
                            BYVAL pa_VEC AS UDT_VEC POINTER , _
                            BYVAL pb_VEC AS UDT_VEC POINTER , _
                            BYVAL pc_VEC AS UDT_VEC POINTER )
                            
                            
                            IF pa_VEC = 00 OR _
                            pb_VEC = 00 OR _
                            pc_VEC = 00 THEN
                            
                            EXIT SUB
                            
                            ELSE
                            
                            @pc_VEC.v(00) = @pa_VEC.v(01) * @pb_VEC.v(02) - @pa_VEC.v(02) * @pb_VEC.v(01)
                            @pc_VEC.v(01) = @pa_VEC.v(02) * @pb_VEC.v(00) - @pa_VEC.v(00) * @pb_VEC.v(02)
                            @pc_VEC.v(02) = @pa_VEC.v(00) * @pb_VEC.v(01) - @pa_VEC.v(01) * @pb_VEC.v(00)
                            
                            END IF
                            
                            END SUB
                            
                            
                            SUB Vector_CRS_NORM ( _
                            BYVAL pa_VEC AS UDT_VEC POINTER , _
                            BYVAL pb_VEC AS UDT_VEC POINTER , _
                            BYVAL pc_VEC AS UDT_VEC POINTER )
                            
                            
                            IF pa_VEC = 00 OR _
                            pb_VEC = 00 OR _
                            pc_VEC = 00 THEN
                            
                            EXIT SUB
                            
                            ELSE
                            
                            Vector_CRS pa_VEC , _
                            pb_VEC , _
                            pc_VEC
                            
                            Vector_NORM pc_VEC
                            
                            END IF
                            
                            END SUB
                            where the vector norm is:

                            Code:
                            SUB Vector_NORM ( _
                            BYVAL pa_VEC AS UDT_VEC POINTER )
                            
                            
                            DIM i_ROW AS LONG
                            
                            DIM NORM_a AS DOUBLE
                            
                            IF pa_VEC = 00 THEN
                            
                            EXIT SUB
                            
                            ELSE
                            
                            NORM_a = Scalar_NORM(pa_VEC)
                            
                            IF NORM_a > 0.00# THEN
                            
                            FOR i_ROW = 00 TO 02
                            
                            @pa_VEC.v(i_ROW) /= NORM_a
                            
                            NEXT i_ROW
                            
                            END IF
                            
                            END IF
                            
                            END SUB

                            Comment


                            • #15
                              Try this with B MAT? I think not. This is compact routine to generate a DCM given any Euler angle sequence. This uses a call by DWORD to select the proper rotation transformation based on the E.ijk permutation.

                              Code:
                              SUB Matrix_T_EAR (w_EAR_ijk AS LONG , _
                              BYVAL pa_EAR AS UDT_VEC POINTER , _
                              BYVAL pM_DCM AS DOUBLE POINTER )
                              
                              ' Routine performs a sequenced Euler angle rotation (EAR)
                              
                              DIM w_EAR_i AS LONG
                              DIM w_EAR_j AS LONG
                              DIM w_EAR_k AS LONG
                              
                              DIM pR_i AS DWORD
                              DIM pR_j AS DWORD
                              DIM pR_k AS DWORD
                              
                              DIM L_INIT AS STATIC LONG
                              
                              IF NOT L_INIT THEN
                              L_INIT = %T
                              REDIM M_i (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                              REDIM M_j (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                              REDIM M_ji(00 TO 02, 00 TO 02) AS STATIC DOUBLE
                              REDIM M_k (00 TO 02, 00 TO 02) AS STATIC DOUBLE
                              END IF
                              
                              
                              IF pa_EAR = 00 OR _
                              pM_DCM = 00 THEN
                              
                              EXIT SUB
                              
                              ELSE
                              
                              ' Test for default condition which is the standard 321 sequence
                              ' Note this is done to speed the most often used 321 sequence
                              
                              IF w_EAR_ijk = 321 THEN
                              
                              ' Creates a DCM by Euler angle rotation 321
                              
                              @pM_DCM[00 OF 02, 00 OF 02] = COS(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                              @pM_DCM[00 OF 02, 01 OF 02] = COS(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                              @pM_DCM[00 OF 02, 02 OF 02] = - SIN(@pa_EAR.v(01))
                              
                              @pM_DCM[01 OF 02, 00 OF 02] = - COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(00)) + SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                              @pM_DCM[01 OF 02, 01 OF 02] = COS(@pa_EAR.v(02)) * COS(@pa_EAR.v(00)) + SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                              @pM_DCM[01 OF 02, 02 OF 02] = SIN(@pa_EAR.v(02)) * COS(@pa_EAR.v(01))
                              
                              @pM_DCM[02 OF 02, 00 OF 02] = SIN(@pa_EAR.v(02)) * SIN(@pa_EAR.v(00)) + COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * COS(@pa_EAR.v(00))
                              @pM_DCM[02 OF 02, 01 OF 02] = - SIN(@pa_EAR.v(02)) * COS(@pa_EAR.v(00)) + COS(@pa_EAR.v(02)) * SIN(@pa_EAR.v(01)) * SIN(@pa_EAR.v(00))
                              @pM_DCM[02 OF 02, 02 OF 02] = COS(@pa_EAR.v(02)) * COS(@pa_EAR.v(01))
                              
                              
                              ELSE
                              
                              ' Parse out the EAR sequence ijk into individual components
                              ' Note this will handle any combination including zeros which are trapped below
                              w_EAR_i = FIX( w_EAR_ijk / 100.0#)
                              w_EAR_j = FIX((w_EAR_ijk - 100 * w_EAR_i) / 10.0#)
                              w_EAR_k = FIX((w_EAR_ijk - 100 * w_EAR_i - 10 * w_EAR_j) / 1.0#)
                              
                              #IF 00
                              MSGBOX "w_EAR i,j,k = " + STR$(w_EAR_i) + " " + _
                              STR$(w_EAR_j) + " " + _
                              STR$(w_EAR_k) , %MB_ICONINFORMATION OR %MB_TASKMODAL , $S_MSGBOX_HDR
                              #ENDIF
                              
                              ' Choose the codepointer for the specific rotation based on the EAR ijk indices
                              pR_i = CHOOSE(w_EAR_i, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                              pR_j = CHOOSE(w_EAR_j, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                              pR_k = CHOOSE(w_EAR_k, CODEPTR(Matrix_R_x), CODEPTR(Matrix_R_y), CODEPTR(Matrix_R_z), ELSE 00)
                              
                              ' Test the function pointer and call the respective rotation matrix by code pointer
                              ' Note the null rotation pointer is trapped by applying the identity matrix
                              IF pR_i = 00 THEN
                              Matrix_IDN VARPTR(M_i(00))
                              ELSE
                              CALL DWORD pR_i USING Matrix_R_n_PATTERN(@pa_EAR.v(00), VARPTR(M_i(00)))
                              END IF
                              
                              IF pR_j = 00 THEN
                              Matrix_IDN VARPTR(M_j(00))
                              ELSE
                              CALL DWORD pR_j USING Matrix_R_n_PATTERN(@pa_EAR.v(01), VARPTR(M_j(00)))
                              END IF
                              
                              IF pR_k = 00 THEN
                              Matrix_IDN VARPTR(M_k(00))
                              ELSE
                              CALL DWORD pR_k USING Matrix_R_n_PATTERN(@pa_EAR.v(02), VARPTR(M_k(00)))
                              END IF
                              
                              ' Create Mji = Mj * Mi
                              Matrix_MLT VARPTR(M_j (00)) , _
                              VARPTR(M_i (00)) , _
                              VARPTR(M_ji(00))
                              
                              ' Create M_DCM = Mk * Mji = Mk * Mj * Mi
                              Matrix_MLT VARPTR(M_k (00)) , _
                              VARPTR(M_ji(00)) , _
                              pM_DCM
                              
                              END IF
                              
                              END IF
                              
                              END SUB
                              
                              #IF %F ' Used for testing
                              
                              REDIM M_EAR(00 TO 02, 00 TO 02) AS DOUBLE
                              REDIM a_EAR(00 TO 02) AS DOUBLE
                              
                              ' Create a DCM Euler matrix based on a ijk rotation
                              
                              a_EAR(00) = 10.0# / 180.0# * Number_PI
                              a_EAR(01) = 20.0# / 180.0# * Number_PI
                              a_EAR(02) = 30.0# / 180.0# * Number_PI
                              
                              Matrix_R_DCM_EUL 123 , _
                              VARPTR(a_EAR(00)) , _
                              VARPTR(M_EAR(00))
                              
                              
                              ? Matrix_GET_S(VARPTR(M_EAR(00)))
                              
                              #ENDIF

                              Comment


                              • #16
                                Thank you Sir Macia , PB MAT method has a very small footprint. Is there a cross product for Matrices using PB MAT?

                                Thank you Sir SchRage, I will test your code out on cross product for matrices

                                Comment


                                • #17
                                  Some confusion here. What I posted was a cross product for vectors. Two vectors cross to another vector. A matrix times a vector is a vector. A dot product of vectors is a scalar. A dot product of matrices is the inner product which is a vector.

                                  I don't think there is an analogy to a cross product of matrices.

                                  Anyway, you're welcome. Glad it worked.

                                  Do you need the full library?

                                  Comment

                                  Working...
                                  X