%%
%% DSP Blockset Code Generation Target File
%% S-Function: sdspmsqr
%% Block name: Matrix Square
%%
%% Copyright 1995-2000 The MathWorks, Inc. 
%%
%% $Revision: 1.3 $
%% $RCSfile: sdspmsqr.tlc,v $
%% $Date: 2000/06/15 20:37:45 $


%implements "sdspmsqr" "C"
%include "dsplib.tlc"


%% Function: Outputs ==========================================================
%%
%function Outputs(block, system) Output
    /* DSP Blockset Matrix Square (%<ParamSettings.FunctionName>) - %<Name> */
    %%
    %assign INPORT  = 0
    %assign OUTPORT = 0
    %assign cplx    = LibBlockInputSignalIsComplex(INPORT)
    %assign Awidth  = LibDataInputPortWidth(INPORT)
    %assign Ac      = SFcnParamSettings.ACOLS
    %assign Ar      = Awidth / Ac
    %%
    %<ErrorIfDiscontiguous(block,INPORT)>\
    %%
    %if Awidth == 1
      %% ================ A is a scalar
      %%
      /* Input is a scalar */
      %if !cplx
         %assign t = LibBlockInputSignal(INPORT, "", "", 0)
         %<LibBlockOutputSignal(OUTPORT,"","",0)> = %<t> * %<t>;
      %else
        %assign tre = LibBlockInputSignal(INPORT, "", "", "re0")
        %assign tim = LibBlockInputSignal(INPORT, "", "", "im0")
        %<LibBlockOutputSignal(OUTPORT,"","",0)> =
                                %<tre> * %<tre> +
                                %<tim> * %<tim>;
      %endif
      %%
    %elseif Ac==1
      %% ================ A is a column vector
      %%
      /* Sum-of-squares */
      %if !cplx
        %% Real sum-of-squares
        {
            %assign y = LibBlockOutputSignal(OUTPORT, "", "", 0)
            %assign u = LibBlockInputSignal(INPORT, "i", "", 0)
            const int_T Arows = %<Ar>; /* Rows in input matrix */
            int_T i;
            %<y> = 0.0;
            for (i=0; i<Arows; i++) {  /* Real sum-of-squares */
                %<y> += %<u> * %<u>;
            }
        }
      %else
        %% Complex sum-of-squares
        {
            %assign y = LibBlockOutputSignal(OUTPORT,"","",0)
            %assign ure = LibBlockInputSignal(INPORT, "i", "", "re0")
            %assign uim = LibBlockInputSignal(INPORT, "i", "", "im0")
            const int_T Arows = %<Ar>; /* Rows in input matrix */
            int_T i;
            %<y> = 0.0;
            for (i=0; i<Arows; i++) {  /* Complex sum-of-squares */
                %<y> += %<ure> * %<ure>
                     +  %<uim> * %<uim>;
            }
        }
      %endif
      %%
    %elseif Ar==1
      %% ================ A is a row vector
      %%
      /* Outer product */
      %if !cplx
        %% Real outer product
        {
            const int_T Acols = %<Ac>; /* Cols in input matrix */
            real_T *y = %<LibBlockOutputSignalAddr(OUTPORT,"","",0)>;
            int_T i;
            for (i=0; i<Acols; i++) {   /* Real outer product */
                const real_T scale = %<LibBlockInputSignal(INPORT, "i", "",0)>;
                int_T j;
                for (j=0; j<Acols; j++) {
                    *y++ = scale * %<LibBlockInputSignal(INPORT, "j", "",0)>;
                }
            }
        }
      %else
        %% Complex outer product
        {
            const int_T Acols = %<Ac>; /* Cols in input matrix */
            creal_T *y = %<LibBlockOutputSignalAddr(OUTPORT,"","",0)>;
            int_T i;
            for (i=0; i<Acols; i++) {
                const creal_T scale = %<LibBlockInputSignal(INPORT,"i","",0)>;
                int_T j;
                for (j=0; j<Acols; j++) {
                    creal_T A2_val = %<LibBlockInputSignal(INPORT,"j","",0)>;
                    y->re   = CMULT_XCONJ_RE(scale, A2_val);
                    y++->im =-CMULT_XCONJ_IM(scale, A2_val);
                }
            }
        }
      %endif
      %%
    %else
      %% ================ A is a full matrix
      %%
      %if !cplx
        %% Real matrix square
        {
            const int_T Arows = %<Ar>; /* Rows in input matrix */
            const int_T Acols = %<Ac>; /* Cols in input matrix */
            real_T *y = %<LibBlockOutputSignalAddr(OUTPORT,"","",0)>;
            real_T *B = %<LibBlockInputSignalAddr(INPORT,"","",0)>;
            int_T   k = Acols;
            while(k-- > 0) {  /* Multiply real A' by A */
                real_T *A1 = %<LibBlockInputSignalAddr(INPORT,"","",0)>;
                int_T i = Acols;
                while(i-- > 0) {
                    real_T *A2 = A1;
                    real_T *B2 = B;
                    real_T  acc = 0.0;
                    int_T   j = Arows;
                    while(j-- > 0) acc += (*A2++) * (*B2++);
                    *y++ = acc;
                    A1 += Arows;
                }
                B += Arows;
            }
        }
      %else
        %% Complex matrix square
        {
            const int_T Arows = %<Ar>;
            const int_T Acols = %<Ac>;
            creal_T *y = %<LibBlockOutputSignalAddr(OUTPORT,"","",0)>;
            creal_T *B = %<LibBlockInputSignalAddr(INPORT,"","",0)>;
            int_T   k = Acols;
            while(k-- > 0) {  /* Multiply complex A' by A */
                creal_T *A1 = %<LibBlockInputSignalAddr(INPORT,"","",0)>;
                int_T i = Acols;
                while(i-- > 0) {
                    creal_T *A2 = A1;
                    creal_T *B2 = B;
                    creal_T  acc = {0.0, 0.0};
                    int_T   j = Arows;
                    while(j-- > 0) {
                        const creal_T A2_val = *A2++;
                        const creal_T B2_val = *B2++;
                        acc.re += CMULT_XCONJ_RE(A2_val, B2_val);
                        acc.im += CMULT_XCONJ_IM(A2_val, B2_val);
                    }
                    *y++ = acc;
                    A1 += Arows;
                }
                B += Arows;
            }
        }
      %endif
      %%
    %endif

%endfunction %% Outputs

%% [EOF] sdspmsqr.tlc